Shahab Karrari

INTEGRATION OF FLYWHEEL ENERGY STORAGE SYSTEMS IN LOW VOLTAGE DISTRIBUTION GRIDS

Shahab Karrari

# INTEGRATION OF FLYWHEEL ENERGY STORAGE SYSTEMS IN LOW VOLTAGE DISTRIBUTION GRIDS

Shahab Karrari

Integration of Flywheel Energy Storage Systems in Low Voltage Distribution Grids

# Integration of Flywheel Energy Storage Systems in Low Voltage Distribution Grids

by Shahab Karrari

Karlsruher Institut für Technologie Institut für Technische Physik

Integration of Flywheel Energy Storage Systems in Low Voltage Distribution Grids

Zur Erlangung des akademischen Grades eines Doktor-Ingenieurs von der KIT-Fakultät für Elektrotechnik und Informationstechnik des Karlsruher Instituts für Technologie (KIT) genehmigte Dissertation

von Shahab Karrari, M.Sc.

Tag der mündlichen Prüfung: 24. Februar 2021 Hauptreferent: Prof. Dr.-Ing. Mathias Noe Korreferent: Prof. Dr.-Ing. Marc Hiller

Bild Umschlagvorderseite: KIT, Institut für Technische Physik (ITEP) Fotografie: KIT, Amadeus Bramsiepe

#### **Impressum**

Karlsruher Institut für Technologie (KIT) KIT Scientific Publishing Straße am Forum 2 D-76131 Karlsruhe

KIT Scientific Publishing is a registered trademark of Karlsruhe Institute of Technology. Reprint using the book cover is not allowed.

www.ksp.kit.edu

*This document – excluding parts marked otherwise, the cover, pictures and graphs – is licensed under a Creative Commons Attribution 4.0 International License (CC BY 4.0): https://creativecommons.org/licenses/by/4.0/deed.en*

*The cover page is licensed under a Creative Commons Attribution-No Derivatives 4.0 International License (CC BY-ND 4.0): https://creativecommons.org/licenses/by-nd/4.0/deed.en*

Print on Demand 2022 – Gedruckt auf FSC-zertifiziertem Papier

ISBN 978-3-7315-1140-3 DOI 10.5445/KSP/1000139457

# **Acknowledgements**

First and foremost, I would like to express my sincere appreciations to my supervisor Prof. Dr.-Ing. Mathias Noe, for his continuous and consistent support throughout my entire PhD, and for his trust in me to pursue my own research ideas.

Next, I would like to thank the DFG research training group "Energy Status Data" for providing the opportunity to carry out my doctoral studies at KIT, and for establishing a platform for interdisciplinary discussions with other researchers. I am very grateful to all the faculty and graduate members of the research training group for the fruitful discussions and constructive feedback during our meetings. In particular, I would like to thank Nicole Ludwig and Michael Vollmer, who have made parts of this work possible through effective collaborations.

I would like to thank Dr. Giovanni De Carne for his support during the last year in reviewing manuscripts and the dissertation, and providing productive feedback for my work. I also thank Dr. Jörn Geisbüsch for proofreading the manuscripts.

The experiments conducted at KIT's EnergyLab 2.0 would have not been possible without the generous help of Mr. Christian Lange and Mr. Frank Gröner. Christian somehow managed to always find a simple and yet effective solution to any hardware issues faced during experiments.

The grid measurements were collected with the cooperation and kind assistance of Mr. Christian Seiler and Mr. Peter Friedmann.

I would like to thank all my colleagues and officemates at ITEP, in particular, Roland Gyuráki, Ruslan Popov, Fabian Schreiner, Nicolò Riva, Philip Kreideweis, Cyra Neugebauer, Tara Benkel, Dustin Kottonau, and Wescley de Sousa. Thank you for being more than just colleagues and hearing me out in times of trouble.

Lastly, I would like to thank my family for their remote but unwavering support during this rather challenging part of my life.

Karlsruhe, December 14, 2020

*Shahab Karrari*

# **Kurzfassung**

Mit dem Ziel, den Stromsektor zu dekarbonisieren und dem Klimawandel zu begegnen, steigt der Anteil erneuerbarer Energieressourcen in den Energiesystemen rund um den Globus kontinuierlich an. Aufgrund des intermittierenden Charakters dieser Ressourcen kann die Aufrechterhaltung des momentanen Gleichgewichts zwischen Erzeugung und Verbrauch und damit der Netzfrequenz ohne angemessene Maßnahmen jedoch eine Herausforderung darstellen. Da erneuerbare Energiequellen mit Umrichterschnittstellen dem System selbst keine Trägheit verleihen, nimmt gleichzeitig die kumulative Systemträgheit ab, was zu schnelleren Änderungen der Netzfrequenz und Bedenken hinsichtlich der Netzstabilität führt.

Ein Schwungrad-Energiespeichersystem (Flywheel Energy Storage System, FESS) kann schnell große Leistungsmengen einspeisen oder aufnehmen, um das Netz nach einer abrupten Änderung der Erzeugung oder des Verbrauchs zu unterstützen. Neben der schnellen Reaktionszeit hat ein FESS den Vorteil einer hohen Leistungsdichte und einer großen Anzahl von Lade- und Entladezyklen ohne Kapazitätsverlust während seiner gesamten Lebensdauer. Diese Eigenschaften machen das FESS zu einem gut geeigneten Kandidaten für die Frequenzstabilisierung des Netzes oder die Glättung kurzfristiger Leistungsschwankungen auf lokaler Ebene.

In dieser Dissertation wird die Netzintegration eines Hochgeschwindigkeits-FESS auf der Niederspannungsebene aus mehreren Perspektiven untersucht. Zunächst wird das Problem der Platzierung und Dimensionierung eines FESS in Niederspannungsverteilnetzen für Leistungsglättungsanwendungen behandelt. Um den am besten geeigneten Standort für ein FESS zu finden, wird eine datengetriebene Methode zur Abschätzung der relativen Spannungsempfindlichkeit vorgestellt, die auf dem Konzept der Transinformation basiert. Der Hauptvorteil der vorgeschlagenen Methode besteht darin, dass sie kein Netzmodell erfordert und nur Messwerte an den interessierenden Punkten verwendet. Messergebnisse aus einem realen Netz in Süddeutschland zeigen, dass mit dem vorgeschlagenen Ansatz die Netzanschlusspunkte mit einer höheren Spannungsempfindlichkeit gegenüber Wirkleistungsänderungen, welche am meisten von einem durch FESS ermöglichten, glatteren Leistungsprofil profitieren können, erfolgreich zugeordnet werden können. Darüber hinaus wird eine neue Methode zur Dimensionierung von Energiespeichersystemen unter Verwendung von Messdaten eingeführt. Der vorgeschlagene Ansatz erkennt wiederkehrende Verbrauchsmuster in aufgezeichneten Leistungsprofilen mit Hilfe des "Motif Discovery"- Algorithmus, die dann zur Dimensionierung verschiedener Speichertechnologien, einschließlich eines FESS, verwendet werden. Anhand von gesammelten Messdaten aus mehreren Niederspannungsnetzen in Deutschland wird gezeigt, dass die Speichersysteme mit den aus den detektierten Mustern abgeleiteten Charakteristika während der gesamten Messperiode effektiv für ihre Anwendungen genutzt werden können.

Als nächstes wurde ein dynamisches Modell eines Hochgeschwindigkeits-FESS entwickelt und mit experimentellen Ergebnissen in mehreren Szenarien, unter Berücksichtigung der Verluste und des Hilfsenergiebedarfs des Systems, validiert. In den untersuchten Szenarien wurde eine maximale Differenz von nur 0,8 % zwischen dem Ladezustand des Modells und dem realen FESS beobachtet, was die Genauigkeit des entwickelten Modells beschreibt.

Nach Festlegung des erforderlichen Aufbaus wurde die Leistungsfähigkeit eines 60 kW Hochgeschwindigkeits-FESS während mehrerer Frequenzabweichungsszenarien mit Hilfe von Power Hardware-in-the-Loop-Tests beurteilt. Die Ergebnisse der PHIL-Tests zeigen, dass das Hochgeschwindigkeits-FESS sehr schnell nach einer plötzlichen Frequenzabweichung reagiert und in knapp 60 ms die erforderliche Leistung erreicht, wobei die neuesten Anforderungen der Anwendungsregeln für die Frequenzunterstützung auf der Niederspannungsebene erfüllt werden.

Um schließlich die Vorteile des schnellen Verhaltens des FESS für Energiesysteme mit geringer Trägheit zu demonstrieren, wurde ein neuartiger adaptiver Trägheits-Emulationsregler für das Hochgeschwindigkeits-FESS eingeführt und seine Leistung in einem Microgrid mit geringer Trägheit durch Simulationen und Experimente validiert. Die Simulationsergebnisse zeigen, dass die Verwendung des FESS mit dem vorgeschlagenen Trägheits-Emulationsregler die maximale Änderungsrate der Frequenz um 28 % und die maximale Frequenzabweichung um 44 % während der Inselbildung des untersuchten Microgrid reduzieren kann und mehrere zuvor vorgestellte adaptive Regelungskonzepte übertrifft. Der vorgeschlagene Regler wurde auch auf einem realen 60 kW FESS mit dem Konzept des Rapid Control Prototyping implementiert, und die Leistungsfähigkeit des FESS mit dem neuen Regelungsentwurf wurde mit Hilfe von PHIL-Tests des FESS validiert. Die PHIL-Ergebnisse, die die allererste experimentelle Validierung der Trägheitsemulation mit einem FESS darstellen, bestätigen die Simulationsergebnisse und zeigen die Vorteile des vorgeschlagenen Reglers.

# **Abstract**

With the aim of decarbonizing the electricity sector and addressing climate change, the share of renewable energy sources in power systems around the globe is consistently increasing. However, due to the intermittent nature of these sources, maintaining the instantaneous balance between the generation and the demand, and therefore, the grid frequency, can be challenging without adequate measures. In addition, since converterinterfaced renewables do not inherently provide inertia to the system, the cumulative system inertia is simultaneously decreasing, resulting in faster changes in the grid frequency and concerns over the grid stability.

A Flywheel Energy Storage System (FESS) can rapidly inject or absorb high amounts of power in order support the grid, following an abrupt change in the generation or in the demand. In addition to the quick response time, a FESS has the advantage of a high power density and a large number of charging and discharging cycles, with no capacity loss throughout its lifetime. These characteristics make the FESS a well-suited candidate for providing frequency support to the grid or smoothing out short-term power variations at a local level.

The work presented in this thesis studies the grid integration of a high-speed FESS in low voltage distribution grids from several perspectives. First, the problem of allocation and sizing of a FESS in low voltage distribution grids for power smoothing applications is addressed. In order to find the most suitable location for a FESS, a data-driven method for estimating the relative voltage sensitivity is introduced based on the concept of mutual information. The main advantage of the proposed method is that it does not require a grid model, and uses only the measurement values at points of interest. Measurement results from several real distribution grids in southern Germany show that the proposed approach can successfully allocated the grid connection points with a higher voltage sensitivity to active power changes, which can benefit the most from a smoother power profile, using a FESS. Moreover, a new method for sizing energy storage systems using historical measurement data is introduced. The proposed approach detects reoccurring consumption patterns in the measured power profiles using the motif discovery algorithm, which are then used for sizing different energy storage technologies, including a FESS. Using collected measurement data from several low voltage distribution grids in southern Germany, it is demonstrated that the storage systems with the characteristics derived from the detected patterns only can be effectively used for their intended applications during the whole measurement period.

Next, a dynamic model of a high-speed FESS has been developed, which incorporates the losses and auxiliary power requirements of the system. The model is validated with experimental results from a real FESS in several scenarios. A maximum difference of only 0.8 % has been observed between the state of charge of the real FESS and its model during the investigated scenarios, reflecting the accuracy of the developed model.

After establishing the required setup, the performance of a 60 kW high-speed FESS during several frequency deviation scenarios has been evaluated by means of Power Hardware-in-the-Loop (PHIL) testing. The PHIL testing results show that the high-speed FESS responds very quickly and reaches the required power in just under 60 ms following a frequency deviation, while complying with the latest grid code requirements for frequency support at the low voltage level in Germany.

Lastly, in order to demonstrate the advantages of the rapid response of the FESS in lowinertia power systems, a novel adaptive inertia emulation controller for the high-speed FESS has been introduced and its performance has been validated in a low-inertia microgrid using simulations and experiments. Simulation results show that the use of the FESS with the proposed inertia emulation controller can reduce the maximum rate of change of frequency by 28 %, and the maximum frequency deviation by 44 % during the islanding of the studied microgrid. In addition, the proposed controller outperforms several previously reported adaptive control designs in terms of limiting the frequency deviations. The proposed controller has also been implemented on a real 60 kW FESS using the concept of rapid control prototyping, and the performance of the FESS with the new control design has been validated by means of PHIL testing. The PHIL results, which are the very first experimental validation of inertia emulation using a FESS, confirm the simulation results, demonstrating the benefits of the proposed controller in low-inertia power systems.

# **Table of Content**





# **1 Introduction, Motivation, and Scope of Work**

In order to comply with the global efforts to address climate change, more and more renewable energy sources are being integrated into power systems around the world. The global installed capacity of renewables has been steadily increasing over the last decades, reaching approximately 200 GW by the end of 2020 [1]. Renewables are projected to become the largest source of energy for electricity generation worldwide by 2025, surpassing coal, and supplying one-third of the global electricity demand. By 2050, over 80 % of the world's electricity could be supplied by renewables [2]. However, maintaining the stability and reliability of the power systems and the power quality can be challenging with such a high share of renewables.

In Germany, the share of renewables in the electricity supply reached approximately 42 % in 2019, exceeding the target of 35 %, initially planned for 2020 according to the renewable energy act [3]. Solar and wind generation units make up the largest share of renewable generation in Germany, with a current total installed capacity of 45 GW and 59 GW, respectively. By 2034, the net installed capacity of wind and solar generation is expected to reach 173 GW, which is approximately twice the current electricity consumption peak in the country. However, due to the intermittent nature of these resources, the instantaneous mismatch between the generation from renewables and demand is expected to grow as high as 74.3 GW [4]. Therefore, balancing the volatile generation and the demand in each and every moment can be challenging without adequate measures. This balance is reflected in the grid frequency, which increases from the nominal value in case of generation exceeding demand, and decreases in case of insufficient generation. For large deviations of frequency above a certain limit, renewables are required to be curtailed, losing useful energy, while for large underfrequency incidents, large consumers are required to be disconnected by the underfrequency load shedding schemes, leading to significant financial losses.

Furthermore, contrary to conventional power plants, converter-interfaced generation units, such as wind and solar, do not inherently provide inertia to the power system. Therefore, with the increasing share of converter-interfaced renewables, the total system inertia is steadily decreasing in power systems around the globe, including in Europe [5]. This leads to a higher rate of change of frequency, which can potentially lead to an increase in the number of frequency violations and activation of protection systems based on a high rate of change of frequency. This issue is expected to become more critical by the decommissioning of nuclear and coal power plants in countries such as Germany in the near future. With the reduced system inertia in the interconnected European power system, a system split, similar to the one occurred on November 4, 2006, can have serious consequences for the European power system [5], [6]. In smaller power systems with less interconnections, such as the power system of the UK and Ireland, the increase in the rate of change of frequency has already caused great concerns, and precautionary measures are being adopted [7].

One of the major solutions is the use of Energy Storage Systems (ESS), which can help balancing the intermittent power of renewables and the demand in a wide range of time scales, depending on the technology. Storage systems can play a fundamental role in integrating renewable energy resources and regulating the grid frequency. The International Energy Agency (IEA) has estimated that in order to support the electricity sector decarburization plans, at least 310 GW of additional energy storage systems are required only in the four major regions of Europe, United States, China, and India [8]. High investment costs and regulatory and market constraints are considered as the main barriers towards a more widespread use of ESS in power systems [9]. For balancing the generation and the demand, ESS are complimented by other technical solutions, including flexible demand and generation, sector coupling (e.g., heat, gas), and increasing the power transfer capabilities with additional power corridors (e.g., HVDC).

In order to support the grid frequency in power systems with a reduced system inertia, an ESS should be able to respond promptly to frequency deviations caused by a sudden change in the renewable generation or in the demand. Therefore, a fast-reacting ESS with a high power density is required, which can rapidly provide or absorb high amounts of power in order to support the grid during disturbances [7], [10]–[12]. While a handful of storage technologies can provide fast frequency support services to the grid, several studies have shown that Flywheel Energy Storage Systems (FESS) are the most costeffective solution for high-power short-term ancillary services such as primary frequency support, due to their low operational and maintenance costs and long lifetime [13], [14]. Supercapacitors [11] and lithium-ion batteries [15] are the main competitors of the FESS for this application. However, contrary to the FESS, they can degrade with time and temperature [16], [17], and often particular charging strategies may be required to be adopted to preserve their lifetime [15]. However, a FESS can quickly provide or absorb high amounts of power with no concern over its lifetime or capacity, and support the grid during sudden variations in generation or in demand. With a high number of cycles in its lifetime, a FESS can also be used for smoothing out short-term power fluctuations of renewables or loads at their connection point to the grid, and avoid fast voltage changes caused by active power changes, and the spread of fast power fluctuations to the rest of system.

The focus of this work is to study the integration of a FESS in low voltage distribution grids, where a significant share of distributed and renewable energy resources is connected. Grid integration of new power system components are commonly investigated using either numerical simulations or experimental tests. With the increasing capabilities and wide-spread availability of real-time simulators, Hardware-in-the-Loop (HIL) and Power Hardware-in-the-loop (PHIL) testing are becoming one of the major tools for grid integration studies of new power system components and validating their performance. By enabling the possibility of verifying the performance of a real hardware in a simulated grid, PHIL testing combines the advantages of experimental validations with the flexibility and of simulations. This allows to test new system components in various grid conditions, including abnormal situations such as extreme frequency deviations, which are not possible or safe to carry out in a field test, due to the safe operational limits of real power systems. Moreover, since the real hardware is tested in PHIL testing, it can reveal existing issues, which are difficult to detect using only numerical simulations [18].

Recently, the PHIL testing of lithium-ion battery energy storage systems [19]–[21], supercapacitors [22], and hybrid battery-supercapacitor storage systems [23] have been reported in literature for different applications, including frequency support. However, to our knowledge, the PHIL testing of a FESS in AC grids, in particular for the application of frequency support has not yet been conducted. There are only limited experiments carried out using a FESS in shipboard DC power systems, where the FESS is supporting pulsed power loads [24], [25].

This work aims to conduct a comprehensive study on the integration of a FESS in low voltage distribution grids and validating its performance by means of PHIL testing. In particular, the following objectives has been defined:


With the aim of addressing the aforementioned objectives, the content of this thesis is organized as follows. In Chapter 2, the FESS is introduced and its main characteristics in comparison to other storage technologies are discussed. The structure of the FESS, its main components, and its major applications in power systems are presented. Here, the latest developments in using a FESS for different applications are reported, including novel control designs and advancements in system components.

Chapter 3 gives a short introduction to the applications of real-time simulations in power systems, including the PHIL testing of distributed energy resources and energy storage systems. The challenges in real-time simulation of low voltage distribution grids, which is an important prerequisite for the PHIL testing of any component in a low voltage grid, are discussed, followed by a summary of the available solutions and techniques to tackle these challenges. Next, a step-by-step guideline for conducting PHIL testing of new components is provided. Here, the recent efforts in PHIL testing of energy storage systems are also presented. Lastly, the 1 MVA PHIL testing facility at KIT's EnergyLab 2.0 and its main components are briefly described.

The problem of allocation and sizing of a FESS in low voltage distribution grids is addressed in Chapter 4. Firstly, the optimal placement of a FESS is investigated by locating the grid connection points with the highest voltage sensitivity to active power changes. After discussing the limitations of the classical voltage sensitivity calculations, a data-driven approach based on the concept of mutual information [26] is introduced. The allocations of a FESS in a real grid in southern Germany is presented as a use case of the proposed approach. The second part of this chapter deals with the problem of sizing of storage systems using historical measurement data. It introduces a new sizing methodology, which uses reoccurring daily consumption patterns detected by the motif discovery algorithm [27], rather than using an arbitrary day, as commonly found in literature. The proposed method is applied for deriving the characteristics of two types of energy storage systems, including a FESS.

Prior to conducting the PHIL tests of any component, having a model of the hardware under test can be significantly beneficial for designing safe and successful experiments. Therefore, a dynamic model of a high-speed FESS is developed in Chapter 5, which consists of the model of each individual component of the FESS, including the permanent magnet synchronous machine, the voltage source converters, and their controllers. The model is simulated in real-time, which allows testing the FESS model in the same grid models that are later used for the PHIL testing of a real FESS.

In Chapter 6, the performance of a 60 kW/3.6 kWh high-speed FESS with a maximum rotational speed of 45,000 rpm is evaluated by means of PHIL testing in several frequency deviations scenarios. The response of the FESS to frequency disturbances,

including the major frequency incident of August 9, 2019 in the power system of the UK, and its compliance with the latest grid code requirements in Germany are investigated. The description of the PHIL setup assembled and used for the testing the FESS is also provided. The results of the PHIL tests in all investigated scenarios are also used to validate the model of the FESS developed in Chapter 5.

As discussed, the prompt response of a FESS can help maintaining the grid frequency, in particular in low-inertia systems, where any imbalance in the generation and the demand can lead to rapid changes in the grid frequency. Therefore, Chapter 7 proposes the use of a FESS for inertia emulation, and a new adaptive controller for this purpose is introduced. The performance of the proposed design is first evaluated by means of simulations in a low-inertia microgrid, and compared with several previously reported adaptive controllers for inertia support. Next, the controller is implemented on the 60 kW highspeed FESS using the concept of rapid control prototyping, and the performance of the FESS with the proposed controller is verified using PHIL testing of the FESS in a low voltage microgrid in different scenarios.

Finally, Chapter 8 summarizes all the findings, and illustrates the potential applications of the proposed methods, and the future possible research in these areas.

# **2 The Role of Flywheel Energy Storage Systems in Power Systems**

In this chapter, the FESS, its characteristics, and its main applications in power systems are described. For each application, including frequency regulation and power smoothing, the current state-of-the-art together with several real-world implementations are presented. Moreover, the structural differences between the main two types of FESS, i.e., the low-speed and the high-speed FESS, are provided. Lastly, a comparison between various technologies used for the main components of a FESS is given.

# **2.1 Introduction to Flywheel Energy Storage Systems (FESS)**

A FESS is a mechanical energy storage system, which stores energy in the form of kinetic energy by accelerating a large rotating mass with a high inertia, i.e., the flywheel. The stored energy can be extracted from the FESS by deaccelerating the high-inertia rotor. While the concept of a FESS is considered among the oldest storage technologies, recent developments in high tensile strength materials and advanced magnetic bearings have made it possible to increase the energy density and efficiency (above 85 % [28]) of this technology, significantly. The configuration of a conventional FESS is illustrated in Figure 2.1. The term FESS refers to the entire system consisting of a flywheel or a highinertia rotor, an electrical machine, back-to-back bi-directional power converters with a common DC link, and auxiliary components required for running the system.

Figure 2.1: Configuration of a conventional Flywheel Energy Storage System (FESS).

The advantages and disadvantages of the FESS in comparison to other storage technologies are listed in Table 2.1. A FESS has several important advantages, including high power density, fast response time and cycling characteristic, high number of cycles, long lifetime, high controllability, insensitivity to ambient temperature, and environmentally benignity. Moreover, the capacity of a FESS remains consistent throughout its entire lifetime, independent of the way it has been charged or discharged and the ambient temperature. In terms of the specific power and energy, the flywheel technology bridges the gap between the supercapacitors and lithium-ion batteries, as shown in Figure 2.2. While a FESS can have the power capabilities of a supercapacitor, it can have a higher specific energy, making longer discharge times possible. Moreover, a FESS has significantly longer lifetime and a greater number of cycles in comparison to lithium-ion batteries.

Table 2.1: Advantages and disadvantages of a FESS in comparison to other storage technologies.


Figure 2.2: Comparison of the specific power and energy of the FESS with supercapacitors and lithiumion (Li-ion) batteries (adopted from [29]).

Limited storage capacity, high self-discharge rate, and high capital costs are considered as the main barriers towards a widespread use of FESS [30]. Higher energy densities and lower self-discharge rates are achieved by the use of composite fiber materials and advanced magnetic bearings to reach higher rotational speeds, but these technologies increase the production costs as well [13]. In addition, the capacity of a FESS can only be increased to a certain extent, due to the limited tensile strength of commercially available rotor materials, which constraints the maximum rotational speed of the rotor (see subsection 2.3.1 for more details). At present, the investment costs of a FESS are relatively high in comparison to several other storage technologies, which could be due to the limited number of flywheel manufacturers worldwide. A list of current manufacturers of a FESS with the characteristics of their product is provided in Table 2.2. With growing interest in the use of a FESS for different applications, the number of manufacturers can increase, which can potentially reduce the capital costs of the FESS. It is estimated that the capital costs of this storage technology would fall by 35 % by 2030 according to a report by the International Renewable Energy Agency (IRENA) [2]. Safety concerns in the case of system failures, are also of crucial importance when using flywheels, but these risks are reduced when light composite fiber materials are used for the rotor [31] (see subsection 2.3.5). While lithium-ion batteries will remain the dominant storage technology for the new storage installations worldwide [32], for applications which require short-term storage of high amounts of power with a high number of cycles, other technologies such as FESS can be more suitable [33].


Table 2.2: A comparison of the characteristics of several commercial FESS [34].

# **2.2 Applications of FESS in Power Systems**

A FESS can be used for a variety of applications in different sectors, depending on its ratings and its type (see section 2.3). Table 2.3 lists the applications of this storage technology in different sectors and areas. In the field of power systems, a FESS is most commonly used for frequency regulation and smoothing out short-term power fluctuations of wind generation units. In weakly interconnected power systems, such as in microgrids, these two applications can be intertwined, as mitigating the power variability of renewables simultaneously leads to a more stable frequency. In addition, a FESS can also be used for improving the system transient stability, voltage regulation, spinning reserve, black start, reliability enhancement, fault ride-through support, and unbalanced load compensation. However, these applications are less common in comparison to the other applications and are less specific to the flywheel technology. The use of a FESS can be particularly beneficial, where a significant number of high-power short-duration cycles ranging from seconds to several minutes are required.

As seen in Table 2.3, FESS are also used in several other sectors, including in the transportation sector [35], for pulsed power applications in industry [36], and in space satellites [37]. These applications are out of the scope of this thesis, and therefore, not discussed further. In the following, the applications of a FESS in power systems are described briefly together with a review of the latest works in each area.


Table 2.3: Applications of FESS in different sectors and fields.

#### **2.2.1 Frequency Regulation**

As discussed, a FESS is a well-suited solution for supporting the grid frequency during imbalances between generation and demand, due to the high power density, quick response time, high reliability, long lifetime, and high number of cycles with no degradation.

A single flywheel unit can be used in microgrids to support the frequency, or tens to hundreds of units can be connected in parallel in a modular configuration to build a utility-scale system. For instance, the company Beacon Power has built and successfully commissioned two 20 MW plants, each with 200 high-speed flywheels in Stephentown, New York, and in Hazel, Pennsylvania. These plants provide fast frequency regulation services to the grid. For the state of New York, the flywheels provide approximately 10% of the whole state frequency regulation requirement. The field data have confirmed the advantages of these systems in terms of costs, emissions, and operational constraints [38].

The combination of the flywheel technology with a battery for providing frequency support has also gained attention over the last few years. Recently, a 9 MW hybrid flywheel-battery energy storage system has been commissioned in the Netherlands to be used for enhanced frequency control, where the flywheel increases the power capabilities of the battery [39]. A smaller flywheel-battery hybrid storage system with the rating of 480 kW was successfully tested in Ireland, where the flywheel responds faster with a higher ramp rate to frequency deviations in comparison to the battery, by design [7]. Similarly, in [40], authors proposed the use of a 50 MW flywheel-based plant for the northern Chilean interconnected system, which can quickly provide the required power for primary frequency control and improving the transient stability of the power system.

Applications of a FESS in a microgrid, in particular for regulating the frequency, have drawn a significant attention over the last 20 years, both in academia and industry. A considerable number of projects have been conducted globally in this area. Since a FESS can simultaneously supply the microgrid with both active and reactive power with a short response time, it can improve the power quality and the dynamic security of the microgrid, significantly. For example, a 500 kW/5 kWh FESS was designed and installed in the Portuguese island of Flores. It was shown that with the help of the FESS, a stable frequency is feasible, even with a high penetration of renewable energy sources [41].

In literature, several new control structures for a FESS for the application of frequency regulation have been proposed. Authors in [42] propose an enhanced frequency control for a low-speed FESS with the aim of reducing the frequency and the DC-link voltage deviations. In [43], a coordinated control strategy between a FESS and a battery energy storage system has been suggested for a more effective droop-based frequency control in microgrids. A fuzzy logic controller for controlling multiple storage systems including a FESS has been introduced in [44], which maintains the state of charge of each storage system in a certain range, while providing enhanced frequency services. Recently, the concept of inertia emulation using a FESS has been introduced [45], and a virtual synchronous generator based on a FESS has been suggested [46]. More information on these works and inertia emulation using a FESS are provided later in Chapter 7.

#### **2.2.2 Power Smoothing**

The advantages of using a FESS is often investigated in the presence of wind generation units, where the FESS can mitigate the high variability of the power generated by wind turbines. This is particularly important for integrating wind generating units into weak or isolated microgrids. With a high number of cycles and no capacity loss, a FESS is a suitable candidate for power smoothing applications, where a continuous operation of the energy storage systems is often required.

Several real-world projects including the FESS in the Glencore Raglan Mine in Quebec has been reported, where the FESS smooths the power produced by a 3 MW wind turbine. Similarly in the Saint Paul Island in Alaska, a FESS is used to reduce wind power fluctuations and the fuel consumption of two diesel generators [47].

In literature, numerous novel control methods have been proposed for the use of a FESS for smoothing the power generated by wind turbines. Recently, an intelligent controller based on the adaptive wavelet fuzzy neural networks has been introduced to mitigate the intermittency in wind generation using a FESS, and improve the transient stability of the power system [48]. The authors in [49] have suggested an energy function-based optimal control approach for reducing the effect of wind power fluctuations on the power system. Here, it has also been shown that the FESS can also improve the low voltage ride-through capability of the wind turbine by operating as a static compensator (STATCOM) during low voltage incidents. In another study, the flux-oriented control and the direct torque control have been compared for a low-speed FESS, and it has been argued that the direct torque control is a better solution for smoothing wind power generation using a FESS [50]. The authors in [51] have presented a control algorithm for a FESS to simultaneously achieve mitigation of power fluctuation and dynamic stability enhancement for an offshore wind farm. Moreover, a supervisory control unit with shorttime wind speed prediction has been suggested to be implemented on a wind farm equipped with a FESS, which has proven to be effective for reducing power fluctuations and the size of the required FESS [52].

By compensating for the short-term volatility of renewable energy resources, a FESS can also reduce the variability of the power output of other distributed energy resources and storage technologies. For instance, a FESS can reduce the power variability of a diesel generator, and therefore, reduce its fuel consumption and emissions, improve its efficiency and lifetime, and decrease the number of required startups and shutdowns, when the generator is used as a backup system [41], [53]. In Nullagine and Marble Bar in western Australia, a 500 kW FESS from ABB in combination with large PV systems have managed to reduce the fuel consumption of the installed diesel generators by 400,000 liters per year and save 1100 tons of greenhouse gas emissions annually [54].

When used together with a battery in a hybrid configuration, a FESS can significantly help preserving the lifetime and capacity of the battery [55], [56]. It has been demonstrated that a battery-flywheel hybrid system outperforms a battery-supercapacitor combination in terms of reducing the peak current drawn from the battery and the overall system efficiency [57]. In Australia, a 2 MW FESS is currently being used to relief a 3 MW battery from wearing faster for mitigating wind power fluctuations [58].

### **2.2.3 Voltage Regulation**

A FESS, like any converter-based system can also provide reactive power for regulating the grid voltage at its connection point and operate similar to a static compensator (STATCOM). However, in low voltage distribution grids with a low X/R ratio, active power provided by the flywheel can be more effective in regulating the grid voltage, in comparison to reactive power (see subsection 4.1.1).

Dynamic voltage compensation on distribution feeders using a FESS has been studied in [59], where it has been shown that a series-connected FESS is more effective for voltage regulation than a FESS connected in parallel to the grid. A series-connected FESS equipped with a matrix converter has been designed and built in [60] for compensating voltage sags, and in [61], a predictive optimal controller has been suggested for such a system. It is also proposed to use a FESS for unbalanced voltage drop compensation and reducing the effect of asymmetric loads on the grid [62].

### **2.2.4 Other Applications in Power Systems**

As listed in Table 2.3, a FESS can be used for several other applications in power system as well. For instance, it has been suggested to use a FESS for enhancing the fault ridethrough capability of an HVDC connection to a wind farm during AC faults [63]. The use of a FESS with a self-organizing fuzzy neural network controller to increase the power transfer capability and transient stability of the power system has been suggested by the authors in [64]. More recently, it has been suggested to use a FESS to limit the initial high power demand of electric vehicle fast chargers from the grid at the moment of connection [65].

# **2.3 Components of FESS**

Flywheel energy storage systems can be grouped into two main categories, i.e., the lowspeed (less than 10<sup>4</sup> rpm) and high-speed (10<sup>4</sup> – 10<sup>5</sup> rpm) FESS, which are very different in their structure and design, and are often used for different applications. A comparison of the structural differences between the two types of the FESS is presented in Table 2.4. Low-speed FESS are the cheaper alternative, which are suitable when a discharge time of several seconds is adequate and a high self-discharge rate is not a concern. High-speed FESS use magnetic bearings, composite fiber rotor, and a vacuum enclosure to reach higher rotational speeds and storage capacities and a significantly reduced self-discharge rate. However, they are significantly more expensive than the low-speed systems. While low-speed FESS are considered a matured technology, high-speed FESS have been only recently commercialized. Figure 2.3 shows the structure of a commercial high-speed FESS.

In the following subsections, important components of a FESS and the different technologies used for each component are discussed.


Table 2.4: Comparison of the main two types of FESS [66], [67].

Figure 2.3: Structure of the container holding the rotor and the electrical machine in a high-speed FESS (Picture courtesy: Stornetic GmbH [68], used with permission).

#### **2.3.1 Rotor**

The energy of a FESS (EFW) depends only on the rotor inertia (J<sup>f</sup> ) and its rotational speed (ωm) according to Eq. (2.1).

$$\mathbf{E\_{FW}} = \frac{1}{2} \mathbf{J\_{f}} \boldsymbol{\omega}\_{\rm m}^{2} . \tag{2.1}$$

As seen in Eq. (2.1), while the capacity of a FESS can also be expanded by increasing the rotor inertia, reaching higher rotational speeds has a greater impact on the storage capacity. The maximum rotational speed is limited by the strength of the rotor material to withstand the centrifugal forces affecting the rotor at high rotational speeds [28]. This is referred to as the material tensile strength. The maximum stress that the rotor material experiences is proportional to the material density and the square of the rotor speed at the rotor's outer surface [69]. The maximum stress has to be below the tensile strengths of the used material by a safe margin. If the material tensile strength is represented by , the energy density ( ), and specific energy () is calculated using [66]

$$\mathbf{e}\_v = \mathbf{K}\boldsymbol{\sigma},\tag{2.2}$$

and

$$e\_m = \frac{\mathcal{K}\sigma}{\rho}.\tag{2.3}$$

In Eq. (2.2) and (2.3), is the density of the rotor material and K is called the shape factor, which depends on the rotor geometry, and it is a measure of effective material utilization [31]. Equation (2.3) indicates that in order to achieve a higher capacity per unit of mass, a material with a high tensile strength and a low density is required. These are the characteristics of composite fiber materials with glass or carbon, as seen in Table 2.5. Therefore, high-speed FESS use composite fiber materials, which make it possible to reach rotational speeds up to 100,000 rpm. Low-speed FESS commonly use variations of steel such as carbon steel [70], which are significantly cheaper than the composite fiber materials. In a low-speed FESS, since the rotational speed of the rotor is limited by the tensile strengths of steel, the capacity is often increased by using larger bulkier rotors (see the two-ton rotor in [71], as an example).

According to Eq. (2.2) and (2.3), the shape factor K also influences the energy per unit of volume and mass. The shape factor depends only on the geometry of the rotor. The rotor geometry can also be optimized to reach higher energy densities, as shown in [72]. However, in practice, solid and hollow cylinders are commonly used for low-speed and high-speed FESS, respectively.


Table 2.5: Characteristics of several materials used for the rotor in a FESS [28], [66].

While the energy of a FESS depends on the rotor material and geometry, its maximum power depends on the ratings of the electrical machine and power converters, which are described in the following subsections.

#### **2.3.2 Electrical Machine**

The high-inertia rotor of a FESS can be coupled to various types of electrical machines, including different designs of synchronous machines, asynchronous machines, and reluctance machines. Each machine type and design has its own advantages and disadvantages. Table 2.6 presents a comparison between the important characteristics of the most commonly used electrical machines in a FESS.

Asynchronous machines are often used in a low-speed FESS. They have a robust construction, high torque, high reliability, and low costs. However, due the presence of copper losses in the rotor, they cannot be used in the vacuum enclosure of a high-speed FESS, where convection cooling is not available. Therefore, a Permanent Magnet Synchronous Machine (PMSM) is often used in a high-speed FESS, which has no field winding and consequently, no copper losses on the rotor. Moreover, PMSM have high specific power, high efficiency, simple and robust structure, and excellent dynamic performance [73]. However, they carry the risk of demagnetization under high temperatures [74]. They also have higher costs and lower tensile strength comparing to asynchronous machines.

Other types and designs of electrical machines such as switched reluctance machines [75], homopolar machines [76], and axial-flux machines [77] have also been used in a FESS. However, they are rarely used in practice in commercial systems.


Table 2.6: Comparison of different types of electrical machines used in FESS [53], [78].

#### **2.3.3 Power Converters**

The power converters control the power exchange between the high-inertia rotor and the grid. A FESS is most commonly connected to the AC grid via back-to-back converters with a common DC link [79], as shown in Figure 2.1. These converters convert the output of the electrical machine with a variable frequency to the fixed frequency required by the grid. For some specific applications of flywheels, such as in wind turbines [80], in uninterruptable power supplies [81], or in electric cranes [82], the flywheel can also be coupled through the DC link of an existing converter using a separate DC-AC converter, as depicted in Figure 2.4. In some cases, an extra bi-directional DC-DC converter may be required to reach the required voltage level in such configurations [83].

Figure 2.4: Flywheel connected to the common DC link of a wind generator for power smoothing applications.

Power converters in a FESS are often conventional Voltage Source Converters (VSC), which are controlled using Pulse-Width Modulation (PWM). Different types of switches (MOSFETs, IGBTs, etc.) can be used, depending on the blocking voltage and current requirements and the switching frequency. Recently, the use of silicon carbide-based MOSFETs has also been investigated for a FESS [84]. Higher switching frequencies in the machine-side converter leads to a reduction in the current and torque ripples of the electrical machine. For the grid-side converter, this leads to decreased harmonics injected to the grid and a smaller grid filter. But this increases the switching losses as well (except when silicon carbide switches are used). Due to the low stator inductance of the electrical machines often used in a high-speed FESS, an LC filter is often required to connect the machine-side converter to the electrical machine, in order to reduce the current ripples and the machine losses [28], [85].

In terms of the converter topology, two-level converters are most commonly used in practice in a FESS. The use of a diode-clamped three-level topology, also referred to as the neutral-point clamped converters, is also suggested for a FESS. This configuration leads to less harmonic distortion, lower machine losses, additional flexibility of the output voltage, higher efficiency, and faster dynamic response, in comparison to a conventional two-level converters, under the same switching frequency [86]. A matrix converter is another alternative suggested to be used in a FESS [60], which does not have a DC-link capacitor in between, and directly converts AC to AC. This reduces the volume and weight of the converters, significantly. However, matrix converters have high harmonic distortions and more complex control and protection systems. The authors in [87] have presented a FESS with a Z-source converter and claimed that this configuration leads to a higher efficiency. Other novel topologies and switching schemes have been suggested to improve the converter efficiencies in a FESS, such as in [88], [89]. However, such systems are rarely used in industry.

### **2.3.4 Bearings**

The flywheel bearings have to support the weight of the rotor, keep it in position, damp out mechanical vibrations (referred to as stiffness), and allow the free rotation of the rotor with minimum losses. Table 2.7 presents the advantages and disadvantages of various bearing technologies used in a FESS.

Traditionally, a low-speed FESS uses mechanical bearings, which have high friction losses and lubrication and maintenance requirements. The friction losses lead to a high self-discharge rate and a low lifetime in such systems. Magnetic bearings are used as the alternative solution, in which the rotor is suspended using a magnetic field, often in a vacuum enclosure, eliminating friction losses and wearing. Therefore, they require much less maintenance, and significantly lower self-discharge rates and higher rotational speeds can be achieved. Magnetic bearings lack friction losses, but other losses resulting from the geometry and variance in a magnetic field, such as eddy current losses, leakage fluxes (stray flux paths), and hysteresis losses are present. Auxiliary mechanical bearings are still required for the initial rotation of the flywheel, or in a case of a magnetic bearings failure or an overload [67].

Magnetic bearings are categorized into active and passive bearings. In active magnetic bearings, the magnetic field is controlled using a coil and a control system with respect to the position and the movement of the rotor. This provides a high stiffness and damping, and results in a highly reliable system with a long lifetime and low maintenance [28]. However, the extra power and the complex control system required for operating such systems are considered as the disadvantages of this type of bearing. They increase the system costs as well.

Passive bearings include the use of permanent magnets. However, permanent magnets should be used with another type of bearing, as they are inherently unstable [53]. The use of High-temperature Superconductors (HTS) in combination with permanent magnets as flywheel bearings leads to less intrinsic losses, compared to any other types of magnetic bearings [90]. Figure 2.5 shows the structure of a flywheel with HTS bearings. The HTS also resists the movement of the rotor due to the pinning effect in HTS materials, which prevents the motion of flux lines [91]. However, HTS bearings have less levitation force, lower stiffness, and lower critical frequencies, and require a cryogenic cooling system. Therefore, hybrid configurations, combining different types of magnetic bearing were proposed. For instance, HTS bearings can be used to either support the weight of rotor or to stabilize the rotor, and active magnetic bearings can be used to improve the stiffness and the damping [92]. Several demonstrators of a FESS with HTS bearings have been built and successfully tested, including in [90], [93], [94]. However, there is no commercial FESS currently available in the market, which uses HTS bearings.


Table 2.7: Comparison of different bearing technologies used in FESS.

Figure 2.5: Structure of a FESS with of High-temperature Superconductors (HTS) bearings.

Bearing losses have to be calculated when designing a FESS. Losses in a FESS also include the losses in the electrical machine, power converters, and AC filters. Moreover, the power required for supplying the active magnetic bearings, vacuum and cooling pumps, control systems, and other auxiliary system required for running the FESS can also be considered as the system losses.

#### **2.3.5 Enclosure**

The major task of the flywheel enclosure is to contain the fragments of the rotor in case of a disastrous system failure at high rotational speeds and to dissipate its energy content. Low-speed FESS require thick heavy enclosures, as a steel rotor can be forcefully dismantled into only few pieces with translational movements [95]. On the contrary, composite rotors used in a high-speed FESS tend to disintegrate into numerous small pieces with rotational movements, resulting in energy dissipation caused by friction. As a results of these differences, it has been reported that for the same energy content of the rotor, the enclosure weight of a high-speed FESS needs to be only half of the rotor weight, while for low-speed FESS, it should be the double [53]. Experimental crash tests are necessary in order to ensure the system safety in case of failures at high speeds.

In a high-speed FESS, the enclosure is also used to provide a vacuum surrounding the rotor, minimizing the air friction and system losses. In low-speed systems, the vacuum can be partial or a low-pressure light gas such as helium is used, since the friction torque is proportional to the density and pressure of the gas surrounding the rotor [53].

## **2.4 Summary**

In this chapter the FESS was introduced along with its main characteristics and applications. As discussed, a FESS has a high power density, quick response time, and high number of cycles with a constant capacity through its lifetime, independent of the way it has been used. Considering these characteristics, the FESS is well suited for applications such as frequency regulation and smoothing out short-term power variations. For each application, the latest developments in academia and industry were reported. Furthermore, the structural differences between the low-speed and the high-speed FESS were presented, where it was shown that how a high-speed FESS can reach higher energy densities and lower losses with the use of light composite fiber materials and advanced magnetic bearings. Therefore, the high-speed FESS is a more promising alternative for power system applications. Lastly, a detailed overview and comparison of available technical solutions together with newly introduced alternatives for each component of the FESS were provided.

# **3 Real-time Simulation and Power Hardware-in-the-Loop Testing for Grid Integration Studies**

In this chapter, real-time simulations and their applications in power systems are introduced. The challenges in real-time simulation of low voltage distribution grids are discussed, and the solutions and the techniques to reduce the model execution time of low voltage grid models are presented. Next, the concept of Power Hardware-in-the-Loop (PHIL) testing is described as a powerful tool for grid integration studies of new grid components, such as energy storage systems. A step-by-step guideline for conducting PHIL testing is presented. The main concerns associated with PHIL testing, including instability issues are also discussed. Lastly, the 1 MVA PHIL testing facility at KIT's EnergyLab 2.0 is described.

# **3.1 Challenges of Real-time Simulations of Low Voltage Distribution Grids**

In real-time simulations, by definition, the status of the simulated model must be updated at every small simulation time step in order for the model to accurately resemble the behavior of its physical counterpart [96]. In the case of power systems, this requires that the complex and large differential equations of the electromagnetic transient model of the corresponding power system to be solved within a predetermined simulation time step. The simulation time step for power system transient simulations can be as low as tens of microseconds, depending on the application and the power system dynamics to be simulated [97]. All calculations have to be finished within the given time step at all times. This is known as the hard real-time constraint. On the contrary, in non-real-time power system simulations, there are no constraints on the time available to solve the model equations. Also, variable time-step solvers can also be used, reducing or increasing the simulation time step, when necessary. This is not an option for real-time simulations, which can only apply fixed-step solvers [96]. If the real-time constraint is violated at any time during the simulation, an overrun has occurred, as shown in Figure 3.1. Overruns can lead to instability and incorrect results, and therefore, should be avoided [98]. Overruns can be resolved by increasing the simulation time step, simplifying the models, or distributing the computational effort on several cores, when possible, without losing the model accuracy or the real-time granularity.

Figure 3.1: The hard real-time constraint and the definition of an overrun.

Real-time simulations can reduce the simulation time of power systems, in particular for large complex grids with a high penetration level of distributed energy resources. It enables modifying the model parameters, while the model is running, without the need for recompiling and rerunning the model. Since real-time simulation models run at the same rate as the actual system, they are also referred to as a *Digital Twin* of the simulated power system [99]. However, the main application of real-time simulations is its capability for real-time data exchange with a physical component. This concept, which is referred to as Hardware-in-the-Loop (HIL) and Power Hardware-in-the-Loop (PHIL) testing is later discussed in detail in this chapter in section 3.2. A comprehensive overview of the possible applications of real-time simulation in power systems is provided in [100].

There are several commercial and non-commercial real-time simulators available, which are tailored for power system simulations. A detailed comparison of the available realtime simulators, in terms of their components, structure, design, and capabilities are given in [97]. These real-time simulation platforms were initially designed for the simulation of large high voltage transmission grids, with the aim of testing protection systems or new components such as static var compensators [100]. Over time, new challenges rose in low voltage distribution grids, with the introduction of converter-interfaced components, such as photovoltaic systems, energy storage systems, and electric vehicle chargers. These new challenges called for innovative solutions, which required new testing methods such as Hardware-in-the-Loop (HIL) and Power Hardware-in-the-Loop (PHIL) testing. This led to a clear interest from both industry and academia to validate these technologies using these methods. For connecting a real hardware to a simulated grid, running the grid model in real-time is mandatory. However, this can be challenging in case of some low voltage distribution grids. For real-time simulation of transmission grids, the propagation delays in long transmission lines are often used for model decoupling and parallelization. However, such methods cannot be applied to low voltage grids, where lines and cables are much shorter [101]. Moreover, many converter-interfaced distributed energy resources are often connected at low voltage grids. These components can increase the model computational time significantly, if modelled in detail. These challenges and the available solutions for tackling them are discussed in detail in the following subsections.

### **3.1.1 Lack of Decoupling Points for Parallel Computing**

Most computers have a multi-core design nowdays. To be able to use this advantage for real-time simulations, the model must have the capability to split its solution matrix and solve it on parallel cores, without losing its accuracy. Real-time simulation software can take advantage of the natural propagation delays in transmission lines for model decoupling, when these delays are greater than the simulation time step. By assuming a simulation time step of 50 µs for generic power system simulations [97], and 300,000 km/s as the wave propagation speed through lines (speed of light), this requires an overhead line to be longer than approximately 15 km in order to have the desired delay for model decoupling [102]. For underground cables, the minimum line length for model decoupling is reduced to 10 km, by assuming a propagation speed of 200,000 km/s for cables. If such relatively long lines are present in the model, the two grid areas on either side of the line can be solved independently on multiple cores using the Bergeron traveling wave line model. However, in low voltage distribution grids, lines and cables are usually much shorter [101] and other solutions have to be found.

Researchers have proposed methods to go around this issue by, for instance, introducing artificial delays into the system with the same length of the simulation time step [103]. However, such delays compromise the model accuracy, significantly [104], [105]. Another solution is to add a spurious shunt capacitor after a transformer [102]. This results in a decoupling possibility at the substation. However, since such a capacitor does not exist in reality, this method is a deliberate change to the system impedance. Therefore, this method also reduces the model accuracy. It can also cause unrealistic high transient voltage values at the transformer [106].

Another important factor to be considered in real-time simulations is the type of modelling technique chosen for formulating the power system equations. The main two well-known approaches are the nodal and the state-space analysis and each power system simulator usually uses one of these techniques [107], [108]. But authors in [109] have come up with an innovative solution that combines the state-space and the nodal analysis for the simulation of electrical systems, which can solve the model equations faster for real-time applications, with the possibility for parallelization. This method is known as State-Space Nodal (SSN). The key idea of the state-space nodal approach is to introduce major nodes in the system of equations and to use these nodes as virtual decoupling points in the circuit. Then each section can be described by local state-space equations, which are linked together through the nodal analysis. As these parts are almost independent, their iterations can be parallelized on multiple cores without any artificial delays [104]. This method does not introduce any impedance change in the model, so the results are highly accurate. However, parallelization using the state-space nodal method is not as effective as other decoupling method, as parts of the calculations cannot be parallelized, and calculations are not fully independent within a single simulation time step [109].

#### **3.1.2 Aggregation of Distributed Energy Resources**

Low voltage distribution grids often have a high penetration level of converter-interfaced elements. These components have higher levels of complexity, in both structure and control, in comparison to the traditional power system elements, such as lines and transformers. Therefore, a detailed simulation of these components requires a high computation effort, in particular when they have a large number of switches at high switching frequencies. The high computation burden can even be beyond the capabilities of CPU-based real-time simulators and require running the models on an FPGA.

While for some studies a detailed converter model for distributed energy resources may not be necessary [110], for applications such as hardware-in-the-loop testing of converter controllers or assessing high-frequency harmonics, a switching or detailed converter model (see subsection 5.2.2) is required. For such simulations, the simulation sampling frequency (one over the simulation time step) should be at least 50 to 100 times larger than the pulse-width modulation frequency [111]. This means that a detailed converter model with a switching frequency of 10 kHz requires time steps of 1 µs or less. In addition, fixed-step non-iterative solvers used in real-time simulation software can hardly handle the instantaneous switching events in the converters [112]. Thus, lower simulation time steps may be required for maintaining the numerical stability. Therefore, CPU-based real-time simulators, which have a minimum simulation time step of few microseconds, may not be capable of simulating detailed models of such converters, and the use of FPGA-based simulators might be mandatory. There are several switching event compensation methods [111], [113] that can increase the required simulation time step to some extent, however, these methods also have limitations. The simulation time step in FPGA-based real-time simulation can be in the nanosecond range. In this case, cosimulation using CPU and FPGA can also be a good alternative, in which the grid model runs on the CPU, while the converter switching or detailed model runs on the FPGA with a sub-microsecond simulation time step [114]. In [101], a large distribution grid and a microgrid were simulated in real-time using this approach. However, like any multi-rate real-time simulation, the interfacing of the simulations on the CPU and the FPGA can be challenging and can cause numerical instability. This problem can be solved by adding snubber RC circuits or other solutions proposed in [115]. For enhancing the FPGA capabilities for real-time simulation of power converters with a high number of switches and high switching frequencies, several solutions have been proposed in the literature. For instance, authors in [115], have managed to decouple the power electronics circuits by using different discretization methods and combining implicit and explicit solvers. Another solution has been proposed in [116], which allows subsets of switches to be treated independently to reduce the model computational time.

However, depending on the intention of the study, the converter models and their controllers can be simplified. For applications such as optimization, scheduling, primary frequency control, voltage stability analysis, and even DC-link voltage stability, the converter average model can be adequate [110]. An example of using average converter models for the application of frequency and voltage support using energy storage systems is later presented in Chapter 5.

There are also methods to decouple the converter model on multiple CPU cores to reduce the model execution time. This can be done using the following solutions:


### **3.1.3 Presence of Asymmetry and Harmonics**

Another characteristic of the low voltage grids is that the current and the voltage asymmetries and harmonics are more dominant at this voltage level, in comparison to higher voltage grids. Asymmetric and nonlinear loads and generation units generate unbalanced distorted currents in the grids, which can significantly affect the voltage waveforms as well, due to the high equivalent grid impedance at this voltage level.

Considering asymmetries leads to more complex models of lines and cables and the necessity of a careful consideration of node grounding, which includes choosing the grounding model and how to parameterize it.

Simulating harmonics can drastically increase the model complexity, while simultaneously demanding a relatively small simulation time step. It is recommended to have a simulation sampling frequency of at least 20 times the maximum harmonic frequency that needs to be accurately represented [111]. At the same time, the simulation time step must be large enough to solve the complex differential equations of the model. This can result in a limited number of components that can be simulated in real-time.

## **3.1.4 A Summary of the Solutions for Real-time Simulation of Low Voltage Distribution Grids**

In the previous subsections, the techniques and solutions that can help to tackle the challenges faced for real-time simulations of low voltage grids were discussed. A summary of these techniques and other methods proposed in literature is presented in Table 3.1. These methods can help to use optimally the computational capability of the real-time simulator, while preserving the model accuracy. For real-time simulation of medium- to large-scale low voltage distribution grids with a high penetration level of distributed resources, several of these techniques have to be combined in order to avoid overruns and achieve an acceptable model accuracy. Reducing the simulation time step of the grid model, can also help with the stability of a PHIL setup, as shown later in subsection 3.2.2.


Table 3.1: Possible techniques and solutions for real-time simulations of low voltage distribution grids with a high share of distributed energy resources.

Due to the aforementioned challenges, low voltage grid models are often simplified to reduce the model execution time. In [118], authors have simulated the CIGRE low voltage benchmark [123] in real-time, only as a single-phase grid, assuming a fully symmetric operation of the grid. In addition, the distributed energy resources have been represented only by simple PQ sources. A 174-node distribution grid model was reduced to an equivalent 7-bus model in [20], in order to simulate the model in real-time with a simulation step of 50 µs. Another example is presented in [106], where the unbalanced coupled PI-lines are replaced with simple RL branches with shunt capacitors for model simplification of the IEEE 123 node feeder. In [101], a large distribution grid with 615 nodes has been simulated in real-time together with a microgrid using the state-space nodal solver, without any major simplifications. Among the available solutions, the use of state-space nodal solver can significantly extend the possibilities of real-time simulation of large low voltage grids, without model simplifications or compromising the model accuracy. Therefore, this method is applied later in Chapters 6 and 7 to simulate the studied low voltage grids and microgrids in real-time.

In addition to real-time simulations of power systems, real-time simulations of power system components, such as the energy storage systems can accelerate their design, analysis and validation [96], [100]. This is discussed later in Chapter 5, where a real-time simulation of a FESS for grid integration studies is provided.

# **3.2 Power Hardware-in-the-Loop (PHIL) Testing of Distributed Energy Resources**

Real-time simulations enable a hardware to interact with the simulation models via I/O modules, which is considered a milestone in model-based design, and testing of new components and controllers [96]. In the 20th century, this was only possible with highly costly analog test benches to mimic the real system. Today, digital real-time simulators do the same task at much lower prices and with a higher accuracy and flexibility [107]. Depending on which part of a system is represented by the real hardware and which part by its real-time model, there are different concepts of using real-time simulations, as shown in Table 3.2.


Table 3.2: Different testing concepts using real-time simulations.

In Controller Hardware-in-the-Loop (CHIL) testing, only low-power signals are exchanged between the real-time simulator and the hardware under test, which is suitable for testing controllers and protection relays. Some examples include CHIL testing of energy management systems [124] and distributed voltage controllers [125]. In CHIL testing, the hardware under test receives signals as it would receive in the field and reacts accordingly. Power Hardware-in-the-Loop (PHIL) testing is an extension to the CHIL testing, in which significant electrical power is exchanged with the hardware in real time using power amplifiers. This is an approximate replication of field testing of power system components, at much lower expenses and risks and with a much higher flexibility and repeatability. Moreover, it allows to test the hardware in extreme conditions such as faults or frequency deviations, which is often not allowed or possible in field tests. PHIL testing can be a powerful tool to investigate the behavior of distributed energy resources and energy storage systems in various grid conditions. It has been shown that PHIL testing of distributed energy resources can show existing issues, that it is difficult to see using pure simulations [18]. This can be due to the difficulty in representing system nonlinearities in the model, such as delays, and component parameter variation. A schematic diagram of PHIL testing for a converter-interfaced distributed generation or an energy storage system is shown in Figure 3.2.

Figure 3.2: Schematic diagram of a typical PHIL setup for distributed energy resources.

There are several important steps that must be taken prior to any PHIL testing. These are listed in Table 3.3 and explained in the following subsections, which can serve as a guideline for carrying out PHIL tests. To summarize, firstly, the hardware under test and the grid under study should be modelled and simulated in real-time, and the models should be validated, when possible. The next step is the evaluation of the stability and accuracy of the PHIL setup, in order to avoid physical damage to the PHIL setup and the hardware under test. Finally, PHIL tests can be executed with adequate safety measures.


Table 3.3: A step-by-step guideline for PHIL testing for grid integration studies [118].

#### **3.2.1 Test Case Development**

Firstly, the hardware under test should be modelled and simulated to understand its behavior prior to testing the physical system. This helps to avoid unexpected behavior of the hardware during the tests and to design suitable and safe test scenarios. The level of the details required for this modelling depends on the intended type of study. Although offline simulation is also beneficial, real-time simulation of the hardware under test has the advantage that it can be directly inserted in the real-time simulation model of the grid, in which the actual hardware is later tested. The developed models can be validated with measurements from the real hardware, when possible.

Next, the grid, in which the hardware is planned to be tested, should be implemented in a real-time simulation environment. The details of the grid model, and the simulation step time should be suitable to represent a real power system for the intended study. As mentioned, it is recommended that the simulation sampling frequency should be at least 20 times greater than the maximum frequency of transients that needs to be represented with an acceptable accuracy, which is also known as model bandwidth [111]. Figure 3.3 shows the minimum simulation time step according to the model bandwidth. As seen, the model bandwidth also determines the minimum amplifier bandwidth, which helps selecting the appropriate power amplifier for the PHIL setup. Moreover, the maximum open loop delay for a stable setup is shown according to the model bandwidth. The stability of the PHIL setup is discussed in the following subsection.

Figure 3.3: Real-time simulation and PHIL testing requirements according to the model bandwidth, defined as the maximum frequency of the harmonics or system dynamics to be represented accurately [126].

Lastly, the scenarios for the PHIL testing should be defined. These scenarios should be first simulated with the developed models of the hardware and the grid, preferably, in a real-time simulation environment. By doing so, one can have a good estimate of the behavior of the hardware under test in each test scenario. The delays and latencies associated with the PHIL setup can also be added to the simulation models to increase the simulation accuracy of the PHIL tests. The role of these delays and latencies on the stability and accuracy of the PHIL setup are discussed next.

#### **3.2.2 Stability and Accuracy Evaluation**

In order to replicate the direct connection of the grid and the hardware under test, an interfacing algorithm is required to exchange signals, as demonstrated in Figure 3.4. In one common approach, known as the ideal transformer method, the voltage values from the real-time simulation are amplified and fed to the hardware under test. The consequent current flowing through the hardware are measured and fed back to the simulation to close the loop. In Figure 3.4, eq is the equivalent grid impedance in the real-time simulation environment, and HuT is the impedance of the hardware under test.

The errors, time delays, and nonlinearities in different components of the PHIL setup can lead to system instability and inaccurate results [98]. For instance, the use of a power amplifier introduces additional control loops, errors, time delays, and dynamics in the system and limits the loop bandwidth [127]. Measurement systems can also cause delays due to filtering, A/D and D/A conversions, and signal transmission. A major factor contributing to the stability and accuracy of the setup is the simulation time step of the real-time simulator, where the methods introduced in Table 3.1 can be used to reduce simulation time step and increase the stability of the setup.

Figure 3.4: (a) Field test setup (ideal case). (b) PHIL setup using the ideal transformer method.

Therefore, the stability and accuracy of the PHIL setup with the hardware under test should be evaluated prior to conducting the experiments. Commonly, the stability of the system is analyzed using stability criteria from linear control theory, such as the Nyquist or Bode stability criterion, or using a transient simulation model of the setup. The accuracy of the system can also be evaluated based on the methods proposed in [128].

Researchers have proposed numerous interfacing algorithms for PHIL testing, which also includes methods that are too complex for implementation, with limited applicability. It is common practice to choose a simple interface algorithm that can guarantee stability with a reasonable accuracy. An ideal interfacing algorithm should be accurate with a low settling time and delay, stable in all testing conditions, with a high bandwidth and a simple implementation. In practice, the ideal transformer method and the damping impedance method are most commonly used [118], [129]. These two widely used methods are explained briefly in the following subsections.

#### **3.2.2.1 The Ideal Transformer Method (ITM)**

The Ideal Transformer Method (ITM) is the most commonly used interfacing algorithm for PHIL testing of linear and nonlinear components [18]. This is due to its simple implementation, high accuracy, and the fact that it requires only current measurements as feedback. In this method, the simulated voltage is amplified and fed to the hardware under test and the current flowing through it is measured and fed back to the simulation. A schematic diagram of the ideal transformer method is illustrated in Figure 3.4 (b). From Figure 3.4 (b), the block diagram of a PHIL setup when using the ideal transformer method can be obtained, which is shown in Figure 3.5.

Figure 3.5: Block diagram of a PHIL setup when using the ideal transformer method.

In Figure 3.5, G<sup>F</sup> (s) and G<sup>B</sup> (s) are the feedforward and feedback transfer functions, respectively. These transfer functions can be easily calculated from the transfer function of each individual component of the PHIL setup. The feedforward transfer function G<sup>F</sup> (s) consists of the delay associated with real-time simulator, i.e., the simulation time step T<sup>s</sup> , delays in signal transmission and D/A conversion, and the transfer function of the power amplifier GPA(s). The sum of all delays in the feedforward transfer function is represented by Td1. Therefore, the feedforward transfer function G<sup>F</sup> (s), is calculated as

$$\mathbf{G\_{F}(s)} = \mathbf{G\_{PA}(s)} e^{-\mathbf{T}\_{\mathrm{d}t}s}.\tag{3.1}$$

The feedback transfer function G<sup>B</sup> (s) on the other hand, also includes the ratio of the impedance of the equivalent grid impedance eq to the impedance of the hardware under test HuT. The delays of the current sensor and measurement, signal transmission, and D/A conversion should also be considered, which are all represented by Td2. The feedback measurements can also include a first-order low-pass filter with a cut-off frequency of . Therefore, the transfer function of the feedback loop is calculated using

$$\mathbf{G}\_{\rm B}(\mathbf{s}) = \frac{\omega\_{\rm f}}{s + \omega\_{\rm f} Z\_{\rm H \rm wT}} \mathbf{e}^{-\mathbf{T}\_{\rm d2} \mathbf{s}}.\tag{3.2}$$

For analyzing the stability of the PHIL setup using the ideal transformer method, the open-loop transfer function GOL(s) can easily be calculated by multiplying G<sup>F</sup> (s) and GB (s), as shown in Eq. (3.3). The total loop delay T<sup>d</sup> is defined as the summation of the feedforward and feedback delays, Td1 and Td2. For simplification purposes, the power amplifier transfer function GPA(s) is also represented by a delay function, which its value is integrated into Td. Therefore,

$$\mathbf{G\_{OL}(s) = G\_F(s)G\_B(s) = \frac{\omega\_\mathbf{f}}{s + \omega\_\mathbf{f}}} \frac{\mathbf{z\_{eq}}}{\mathbf{z\_{HuT}}} \mathbf{G\_{PA}(s)e^{-\mathbf{T\_{dz}}s}e^{-\mathbf{T\_{dz}}s} = \frac{\omega\_\mathbf{f}}{s + \omega\_\mathbf{f}} \frac{\mathbf{z\_{eq}}}{\mathbf{z\_{HuT}}} e^{-\mathbf{T\_{dz}}s}.\tag{3.3}$$

A system with the open-loop transfer function of Eq. (3.3) is a conditionally stable system, meaning it is stable only for specific gain values. This can be shown by the Nyquist diagram of the open-loop transfer function GOL(s), as illustrated in Figure 3.6. Using the Nyquist stability criteria has the advantage that the approximation of the time delay, using the Padé approximation [18], is not required. According to the Nyquist stability criteria, a system is unstable, if it encircles the point [-1,0], while the frequency increases on the Nyquist diagram [130]. As seen in Figure 3.6, with impedance ratio below one, the system becomes unstable. Therefore, the ideal transformer method can become unstable, when the ratio of impedance of the hardware under test to the equivalent grid impedance is less than one, i.e., HuT⁄eq < 1.

The parameters used for the stability analysis in this chapter are listed in Table 3.4.

Figure 3.6: The influence of the system impedance ratio on the PHIL setup stability using the ideal transformer method.


Table 3.4: The parameters used for the PHIL stability analysis.

For low values of HuT⁄eq the system can become stable by having a low-pass filter with an adequately low cut-off frequency of the filter [131], [132]. An example is illustrated in Figure 3.7 for a HuT⁄eq of 0.1. As shown, a filter cut-off frequency lower than 2.5 kHz ensures system stability. However, being close to the point [-1, 0] indicates an oscillatory response [130], therefore further reduction of the cut-off frequency may be necessary. This can also be shown by calculating the gain and phase margin of the transfer function using the Bode stability criteria. As shown in Figure 3.8, a positive gain and phase margin is only achieved for a cut-off frequency of below 2.5 kHz. Similarly, since small phase and gain margins indicate an oscillatory closed-loop response, the cutoff frequency should be lower than 2.5 kHz.

Figure 3.7: The influence of the filter cut-off frequency on the PHIL setup stability using the ideal transformer method for a HuT ⁄ eq ratio of 0.1 using the Nyquist criteria.

Figure 3.8: The influence of filter cut-off frequency on the power hardware-in-the-loop setup stability using the ideal transformer method for a HuT ⁄ eq ratio of 0.1 using the Bode criteria.

However, using filters can lead to phase shifts and magnitude attenuations, which reduce the accuracy of the PHIL tests, and can lead to an unacceptable reactive power transfer [18]. Other methods to improve the stability of the ideal transformer method include adding additional components such as inductors [133], using the current-type ideal transformer method [129], interface compensation methods [113], and multi-rate interfacing [134]. Adding inductive components can improve the system stability but increases the hardware impedance, which obviously leads to a reduced accuracy and increased power losses, which is not preferred, in particular for in high-power applications.

The total loop delay T<sup>d</sup> also has a significant effect on the system stability. The most dominant component of the loop delay is the simulation time step of the real-time simulation [135]. As discussed in section 3.1, the simulation time step is constrained by the model complexity and the required model bandwidth that needs to be simulated in real time. Consider the same impedance ratio of HuT⁄eq ratio of 0.1 and a cut-off frequency of 1 kHz, which was shown to be a stable system for a simulation time step of 10 µs. Figure 3.9 depicts the phase and gain margins of the closed-loop transfer function as a function of the simulation time step, ranging from 0 to 100 µs. As seen, the system becomes unstable for simulation time steps of more than approximately 42 µs. Smaller simulation time steps lead to higher stability margins. Since using filters can reduce the model accuracy as previously discussed, reducing the simulation time step should be a priority for improving overall system stability.

Figure 3.9: The influence of simulation time step on the power hardware-in-the-loop setup stability using the ideal transformer method for a HuT ⁄ eq ratio of 0.1 using the Bode criteria.

#### **3.2.2.2 Damping Impedance Method (DIM)**

When the ideal transformer method cannot provide the adequate stability, other interfacing algorithms should be used. The Damping Impedance method (DIM) is another commonly used interfacing algorithm for PHIL testing, which can provide a higher stability margin in comparison to the ideal transformer method [129], [136]. A schematic of a PHIL setup using the damping impedance method is shown in Figure 3.10.

Figure 3.10: Schematic diagram of the Damping Impedance Method (DIM).

As shown in Figure 3.10, a linking impedance <sup>12</sup> is added to both software and hardware side of the PHIL setup, while a damping impedance ∗ is added only to the real-time simulation. Figure 3.11 shows the block diagram of a PHIL setup using the damping impedance method. In Figure 3.10 and Figure 3.11, FF is the voltage feedforward transfer function, which includes the delay associated with signal transmission Td1 and the power amplifier GPA(s), while FB1 and FB2 are the transfer function of the current and voltage feedbacks, respectively.

Figure 3.11: Block diagram of PHIL setup when using the damping impedance method.

The overall feedforward transfer function G<sup>F</sup> (s) using the damping impedance method can be easily calculated by considering the voltage divider shown on the right-hand side of the Figure 3.10. Therefore,

$$\mathbf{G\_{F}(s)} = \frac{\mathbf{v\_{2}}}{\mathbf{v\_{1}}} = \frac{\mathbf{z\_{HuT}}}{\mathbf{z\_{12}} + \mathbf{z\_{HuT}}} \mathbf{G\_{FF}(s)} = \frac{\mathbf{z\_{HuT}}}{\mathbf{z\_{12}} + \mathbf{z\_{HuT}}} \mathbf{G\_{PA}(s)} e^{-\mathbf{T\_{d1}s}}.\tag{3.4}$$

Calculating the overall feedback transfer function G<sup>B</sup> (s) is a bit more complex, as both the voltage and current measurements must be considered. As shown in [129], the overall feedback transfer function G<sup>B</sup> (s) is calculated using

$$\mathrm{G\_B(s)} = \frac{\omega\_f}{s + \omega\_f} \frac{\left(\mathrm{G\_{FB1}\frac{2^\*Z\_{\mathrm{eq}}}{Z\_{\mathrm{H4T}}}\mathrm{G\_{FB2}Z\_{\mathrm{eq}}}}{\left(Z\_{\mathrm{eq}} + Z\_{12} + Z^\*\right)}\right)}{\left(Z\_{\mathrm{eq}} + Z\_{12} + Z^\*\right)}.\tag{3.5}$$

By aggregating all system delays into one variable T<sup>d</sup> and assuming that both voltage and current feedback filters are identical with the same low-pass filter, the open-loop transfer function GOL using the damping impedance method is calculated by

$$\mathbf{G}\_{\rm OL}(\mathbf{s}) = \mathbf{G}\_{\rm F}(\mathbf{s})\mathbf{G}\_{\rm B}(\mathbf{s}) = \frac{\omega\_f}{s + \omega\_f} \frac{\mathbf{z}\_{\rm eq}(\mathbf{z}^\* - \mathbf{z}\_{\rm HT}\mathbf{r})}{(\mathbf{z}\_{12} + \mathbf{z}\_{\rm Ht}\mathbf{r})(\mathbf{z}\_{\rm eq} + \mathbf{z}\_{12} + \mathbf{z}^\*)} \mathbf{G}\_{\rm PA}(\mathbf{s})\mathbf{e}^{-\mathbf{T}\_{\rm dS}}.\tag{3.6}$$

From Eq. (3.6), it is clear that when <sup>∗</sup> = HuT, the magnitude of the open-loop transfer function becomes zero, which indicates a stable system, regardless of all the other variables such as the impedance ratio. However, this improved stability is achieved only when the damping impedance ∗ is equal to the equivalent impedance of the hardware under test HuT. Therefore, the damping impedance method requires a good estimate of the impedance of hardware under test. This can be challenging in practice, in particular for active nonlinear components as hardware under test. In [129], [137], authors propose the use of online impedance identification techniques to match the simulated damping impedance with the impedance of the hardware under test. However, these methods are only applied to complex loads as the hardware under test and their application of power electronic-based devices is not yet investigated. They also add additional hardware and software complications to the PHIL setup. The damping impedance method is difficult to implement in comparison to the ideal transformer method and requires feedback from voltage measurements, which can introduce noise into the real-time simulation. Moreover, adding the linking impedance on the hardware side leads to power losses, which can be problematic for high-power applications. Therefore, this method should be used only if the ideal transformer method cannot provide a stable setup.

#### **3.2.2.3 A Comparison of the Common Interfacing Algorithms**

As discussed previously, there are several interfacing algorithms for PHIL testing proposed in literature, many of which are not discussed in the chapter. Each method can demonstrate improvements in terms of stability and accuracy in particular conditions, while often adding more complications to the setup. Table 3.5 provides a short comparison of the most commonly used interfacing algorithms.


Table 3.5: A comparison of the most commonly used interfacing algorithms for PHIL testing [98], [129], [138].

### **3.2.3 Execution of the PHIL Testing**

After developing the real-time simulation models and the test scenarios and a careful evaluation of the stability and accuracy of the PHIL setup, the hardware under test can be connected to the PHIL setup, using the chosen interfacing algorithm.

Even though a theoretical stability analysis can provide a valuable insight on the stability margin of a PHIL setup, care must be taken to avoid damage to the PHIL facility or the hardware under test, when unexpected instabilities occur [98], [113]. Errors in modelling and parameterization of the PHIL components can lead to an inaccurate stability analysis. For instance, it has been shown in [139] that the total loop delay can vary in time in a PHIL setup, as a result of the interaction between the discrete components of the setup. Therefore, hardware and software protective measures that can shut down the entire system, in case of instabilities or unacceptable high voltages and currents are mandatory.

For executing the PHIL tests, first the grid model should be running in real time. It is then recommended to close the loop with the measurement feedback when there is no or negligible power transfer between the hardware under test and power amplifiers. This can avoid hardware damages in case of incorrect settings or instability. When the loop is closed and the system shows a stable condition, the previously defined scenarios can be executed.

# **3.2.4 State-of-the-art in PHIL Testing of Energy Storage Systems**

PHIL testing has been used to investigate the behavior of the different types of distributed energy resources under realistic conditions, including PV systems [98], [140]–[142], wind turbines [143], and energy storage systems. A summary of the literature review on the PHIL testing of energy storage systems is provided in Table 3.6.

In [20], PHIL testing of an inverter for a battery energy storage systems has been conducted in order to investigate the storage functionalities such as peak shaving and volt-watt control in a distribution grid. An adjustable DC source has been used in this work to emulate the battery. The applications of primary frequency control using a battery energy storage system has been tested in [19] using PHIL testing. It has been shown that the battery has a short response time of around 80 ms. A supercapacitor energy storage system has also been tested in a PHIL environment for frequency support in isolated areas [22]. Batteries and supercapacitors are combined to form a hybrid storage system and tested using PHIL testing in [23], where the supercapacitor is used to smoothen the high power gradients, whereas the battery provides power in the long term.

In case of PHIL testing of a FESS, limited experiments have been done with the focus on shipboard power systems for supporting pulsed loads in DC grids in [24], [25]. To our knowledge, PHIL testing of a FESS in low voltage distribution grids and testing grid functionality such as frequency support using PHIL testing has not yet been conducted.


Table 3.6: Literature review on PHIL testing of energy storage systems.

# **3.3 Description of the 1-MVA PHIL Infrastructure at KIT**

As of December 2019, a 1-MVA PHIL testing facility has made operational by KIT-ITEP, as part of the KIT's EnergyLab 2.0.

An overview of this PHIL setup is illustrated in Figure 3.12. The core of this setup is an Opal-RT 5700 real-time digital simulator with 8 Intel Xeon processing cores and a Xilinx Virtex-7 FPGA on a VC707 board. In addition to the multiple analog and digital I/O boards on the simulator, the system also has 16 high-speed fiber-optic Small Form-factor Pluggable (SFP) sockets (up to 5GBps), which are used for connection to the amplifiers and to multiple OP4520 Kintex7 FPGA & I/O expansion units using Opal-RT's Multi-System Expansion (MUSE) link.

The software RT-LAB is used to compile the models in C programming language and load and execute the models on the real-time target. RT-LAB is also used for monitoring the real-time execution of the models, checking for overruns, assigning I/Os to model inputs and outputs, and changing model parameters, while the model is running on the target.

Figure 3.12: Overview of the 1 MVA PHIL setup at KIT's EnergyLab 2.0.

The voltage set points from the simulated grid are sent digitally to five 200 kVA switched mode amplifiers using the Xilinx's Aurora protocol. The use of this fast digital communication reduces the loop delay caused by signal transmission, by eliminating the D/A conversion unit and their corresponding anti-aliasing filters, and also the problems associated with noise and the signal grounding. The amplifiers are five 200 kVA COMPISO System Unit (CSU) GAMP6 from Egston Power. The nominal characteristics of each amplifier of the PHIL setup is presented in Table 3.7, and the inside view of each amplifier cabinet is shown in Figure 3.13. Each amplifier consists of 6 COMPISO Digital Amplifier (CDA), which can be connected and controlled in various configurations. To simulate a 4-wire three-phase power system, such as the low voltage grids in Europe, they are operated in 3-phase plus neutral mode, in which each phase composes of one CDA, and the neural has 3 CDAs for a high-power neutral connection.

The five Egston amplifier (CSUs), shown in Figure 3.14, can be operated in parallel or in series with each other to reach higher voltage and current ratings with a maximum power rating of 1 MVA. In this case, the set points are first sent an Egston control unit and from there to each individual amplifier.


Table 3.7: The characteristics of each Egston amplifier at KIT-ITEP.

Figure 3.13: Inside of a 200 kVA COMPISO System Unit (CSU) power amplifier from Egston Power [126]. Cabinet 1 shows the grid connection and the isolation transformer, while the grid-side converter and its filter is in cabinet 2. The 6 COMPISO Digital Amplifier (CDA) are in cabinet 3, along with the Egston controller. Cabinet 4 shows the output contactors for the hardware under test.

Figure 3.14: Five 200 kVA GAMP6 Egston switched-mode power amplifiers (1 MVA in total).

For closing the loop for PHIL testing, the current and voltage measurements of the hardware under test are sent to two OP4520 Kintex7 I/O expansion units, which are equipped with an analog input card with 16 channels, 16 bits each, with 2.5 µs conversion time. The measurement values are also sent using high-speed fiber optics communication to the real-time simulator. Currents are measured using three categories of current transformers from the company Danisense, which can measure currents up to 640 A, 2 kA and 5 kA, respectively.

An emergency shutdown system has also been incorporated in the 1 MVA PHIL setup, which opens the amplifiers output breakers and de-energize the system, when necessary. In addition to the hardware protection schemes, software protection algorithms are also implemented on the real-time simulation to open the loop, send trip signals to the amplifier and hardware under test and force the set points to zero, when instable conditions or extremely high currents or voltages are observed.

# **3.4 Summary**

This chapter introduced real-time simulations and its applications in power systems, including PHIL testing. For conducting PHIL tests, a real-time simulation of the grid under study is required. Therefore, a comprehensive overview of the challenges of realtime simulation of low voltage distribution grids, and the available techniques and solutions for tackling these challenges was provided. As discussed, methods that do not comprise the model accuracy such as using faster system solvers should be preferred, and often several methods should be combined to make the real-time simulation of complex grid models feasible.

Moreover, a guideline for conducing PHIL testing was presented. As discussed, having a model of the hardware under test and selecting an appropriate interfacing algorithm are necessary steps for designing safe and fruitful experiments. Among the available interfacing algorithms, the ideal transformer method should be preferred, whenever it can provide a stable setup. This is due to its simple implementation, high accuracy, and the fact that it does not require extra hardware components to be added, which can be a decisive factor for high-power PHIL setups. Moreover, it was shown that a reduction in the simulation step time of simulated grid can reduce the overall loop delay of the setup, and increase the stability margin.

Lastly, a literature review on PHIL testing of energy storage systems revealed that while such tests have been carried out for several storage technologies, a PHIL testing of a FESS in AC grids for applications such as frequency support has not yet been conducted.

# **4 Data-driven Allocation and Sizing of Energy Storage Systems**

Among the major questions concerning the use of an ESS in distribution grids is how to determine its optimal location and characteristics, such as its capacity. It is desired to have the ESS in a location of the grid, where it can be the most beneficial to its stakeholders. After the location of the ESS is determined, it is also necessary to ensure that the installed ESS has the required characteristics to fulfil its main purpose of installation at that particular location. The main characteristics of an ESS include its rated capacity, rated power, and maximum ramp rate. Extensive literature has been published on the topic of optimal allocation and sizing of an ESS [148]–[150]. However, the proposed solutions can vary significantly depending on the storage technology, the intended application of the ESS, and the available grid and measurement data. Therefore, they are not always applicable to other use cases and applications. With the increasing amount of measurement data from low voltage distribution grids, data-driven techniques can provide new insights and solutions to the problem of allocation and sizing of an ESS.

In this chapter, two novel data-driven methods are introduced for the allocation and sizing of storage systems in low voltage distribution grids. Firstly, a new method for estimating the relative voltage sensitivity is introduced based on the concept of Mutual Information (MI) [26] from information theory. The method has been applied to find the most suitable location for a FESS for power smoothing applications in a large interconnected grid in southern Germany, with more than 1200 possible installation candidates. Next, a new method for sizing ESS based on historical measurement data is introduced, that uses reoccurring daily consumption patterns for the sizing, which are detected using the motif discovery algorithm [27]. The detected consumption patterns have been used for finding the characteristics of two types of storage technologies to form a hybrid ESS: A FESS for the application of power smoothing, and a high energy density ESS such as a Li-ion battery, for the application of peak shaving and load balancing. The sizing methodology has been applied for sizing these two types of storage technologies at four different low voltage distribution grids with different levels of PV penetration, and using data with different resolutions.

## **4.1 A Data-driven Method for Allocation of Flywheel Energy Storage Systems<sup>1</sup>**

The optimal location for an ESS in a grid depends strongly on the main objective for its installation and the characteristics of the storage technology. An ESS can have various applications and functionalities in power systems. Therefore, there are numerous techniques proposed in literature for finding the optimal location of different types of ESS in power systems for different applications, including in [151]–[154]. Many of such methods are summarized and compared in [148]–[150]. However, the proposed solutions often focus on optimal allocation of battery energy storage systems for applications such as peak shaving, energy arbitrage, and voltage regulation, and they cannot necessarily be applied for other applications and storage technologies. The work focuses on allocating a high-speed FESS for the application of power smoothing in distribution grids.

A FESS has the ability to quickly and continuously smooth out short-term power fluctuations, as demonstrated in [155], [156], and [157], with no concern over its lifetime or capacity. Short-term power fluctuations in power systems can be caused by variation in loads, for example, industrial loads, or by intermittent generation from renewables, for instance by PV generation units. Power fluctuations at a relatively low amplitude are generally not a concern for the grid operator, but large power fluctuations could lead to significant voltage fluctuations in weak parts of the grid [158]. Fast voltage changes can be a power quality problem according to the European grid code DIN EN 50160 [159]. This is particularly the case in low voltage distribution grids with a low X/R ratio. As shown later in subsection 4.1.1, in such grids, the active power has as much as or even greater impact on the voltage, as the reactive power. In grids with a low X/R ratio, reactive power compensation can be inadequate for voltage regulation, unless an oversized inverter is used [160], while active power compensation can significantly help to mitigate voltage deviations [161], [162]. A FESS can help smoothen the active power variations, which can potentially reduce the voltage fluctuations caused by the active power changes [163].

A FESS can be most useful at the grid connection points, where the voltage is more sensitive to active power changes. This can be quantified using the voltage sensitivity formulation, which is described in detail in the next subsection. The higher is the voltage sensitivity to active power variations, the higher is the impact of active power smoothing using a FESS on regulating voltage. The use of voltage sensitivity has been previously

<sup>1</sup> This section presents a revised version of the following publication: S. Karrari, M. Vollmer, G. De Carne, Klemens Böhm, Jörn Geisbüsch, "A Data-driven Approach for Estimating Relative Voltage Sensitivity," in the 2020 IEEE Power & Energy Society General Meeting, Montreal, 2020.

proposed for the optimal allocation of distributed energy resources to improve the system reliability indices [164], to increase the voltage stability margins [165], and to avoid over- and under-voltage incidents [166]. For the purpose of post-fault voltage recovery, a battery placement method based on voltage sensitivity has shown to be the most effective method in comparison to several other allocation strategies [167], [168]. Voltage sensitivity coefficients have also been used for designing advanced controllers for PV systems in order to improve the PV hosting capacity [169]. However, estimating voltage sensitivity in the aforementioned works requires a validated grid model, which is not always available, in particular for low voltage distribution grids.

This work proposes a new data-driven method for estimating the relative voltage sensitivity to active power changes, which is applied for finding the most suitable location for a high-speed FESS for power smoothing applications. The proposed method requires only measurement values at the points of interest, and does not require a grid model, which is a significant advantage. The suggested approach can also be used for other storage technologies and investigations, that require identification of the locations in a grid with a high voltage sensitivity to active or reactive power changes.

In the following subsections, after a short introduction to voltage sensitivity, the shortcomings of the classical method to calculate analytically the voltage sensitivity coefficients based on steady-state power flow equations are highlighted. Then, a novel data-driven approach has been proposed for estimating the relative voltage sensitivity. The proposed approach has been applied for finding the optimal location for a FESS among more than 1200 10/0.4 kV substations in a region in southern Germany, using the measurement data collected at the top candidates for the installation of the FESS.

#### **4.1.1 Introduction to Voltage Sensitivity**

Figure 4.1 shows a simplified power system with an equivalent impedance, consisting of an equivalent resistance (Req) and an equivalent reactance (Xeq) between the point of common coupling (Bus r) and the main supply (Bus s). All values are assumed to be in per unit.

Figure 4.1: The single-line diagram of a simplified power system for calculating the voltage sensitivity.

The apparent power drawn from Bus r can be calculated as

$$S\_{\mathbf{r}} = P\_{\mathbf{r}} + jQ\_{\mathbf{r}} = V\_{\mathbf{r}}I\_{r}^{\*} = V\_{\mathbf{r}} \left( \frac{\nu\_{\text{s}}e^{j\beta} - \nu\_{\mathbf{r}}}{\mathbf{R\_{eq}} + j\mathbf{X\_{eq}}} \right)^{\*}. \tag{4.1}$$

By multiplying both the nominator and the denominator of the fraction on the right hand side of Eq. (4.1) by (Req − Xeq), and separating the real and imaginary components, Eq. (4.1) can be divided into [170]

$$P\_{\rm r} = V\_{\rm r} \left( \frac{\mathcal{R}\_{\rm eq} (V\_{\rm s} \cos(\delta) - V\_{\rm r}) + \mathcal{X}\_{\rm eq} V\_{\rm s} \sin(\delta)}{\mathcal{R}\_{\rm eq}^2 + \mathcal{X}\_{\rm eq}^2} \right) \tag{4.2}$$
 and

and

$$Q\_r = V\_r \left(\frac{X\_{\rm eq} \langle V\_{\rm s} \cos \langle \delta \rangle - V\_r \rangle - R\_{\rm eq} V\_{\rm s} \sin \langle \delta \rangle}{R\_{\rm eq}^2 + X\_{\rm eq}^2} \right). \tag{4.3}$$

It can be seen from Eq. (4.2) and (4.3) that the relationship between the voltage and the active and reactive power is nonlinear. Therefore, linearization around a single operating point is required to calculate the voltage sensitivity to active and reactive power changes. To calculate the voltage sensitivity to a small disturbance of <sup>r</sup> and <sup>r</sup> around the operating point (r0,r0, r0), Eq. (4.2) and (4.3) are linearized around this point. For doing so, <sup>r</sup> is replaced by r0 + , <sup>r</sup> by r0 + , and <sup>r</sup> by r0 + , and is eliminated, with the trigonometric method described in [170]. The second-order terms such as <sup>r</sup> 2 and <sup>r</sup> 2 are neglected. Therefore, for instance, <sup>r</sup> 2 is replaced by r0 <sup>2</sup> + 2r0 . After these mathematical manipulations, r r and <sup>r</sup> r can be easily calculated, which are assumed to be the approximations of <sup>r</sup> r and <sup>r</sup> r for small changes in the power around this operating point. Hence,

$$\frac{\partial V\_{\rm r}}{\partial P\_{\rm r}} \approx \frac{\Delta V\_{\rm r}}{\Delta P\_{\rm r}} = \frac{\left\{\mathbf{R}\_{\rm eq}^2 + \mathbf{X}\_{\rm eq}^2\right\} P\_{\rm r0} + \mathbf{R}\_{\rm eq} V\_{\rm r0}^2}{V\_{\rm r0} \left(V\_{\rm s}^2 - 2V\_{\rm r0}^2 - 2\left\{\mathbf{R}\_{\rm eq} P\_{\rm r0} + \mathbf{X}\_{\rm eq} Q\_{\rm r0}\right\}\right)}\tag{4.4}$$

and

$$\frac{\partial V\_{\rm r}}{\partial Q\_{\rm r}} \approx \frac{\Delta V\_{\rm r}}{\Delta Q\_{\rm r}} = \frac{\left\{ \mathbf{R}\_{\rm eq}^2 + \mathbf{X}\_{\rm eq}^2 \right\} \mathbf{P}\_{\rm r0} + \mathbf{X}\_{\rm eq} V\_{\rm r0}^2}{V\_{\rm r0} \left( V\_{\rm s}^2 - 2 V\_{\rm r0}^2 - 2 \left\{ \mathbf{R}\_{\rm eq} P\_{\rm r0} + \mathbf{X}\_{\rm eq} Q\_{\rm r0} \right\} \right)}. \tag{4.5}$$

The outputs of Eq. (4.4) and (4.5) are referred to as voltage sensitivity coefficients to active and reactive power changes, respectively, around the operating point (r0,r0, r0). By assuming <sup>s</sup> and r0 in these equations to be approximately 1 per unit, it results in

$$\frac{\partial V\_{\rm r}}{\partial P\_{\rm r}} \approx -\frac{\{\mathbf{X}\_{\rm eq}^{2} + \mathbf{R}\_{\rm eq}^{2}\}P\_{\rm r0} + \mathbf{R}\_{\rm eq}}{1 + 2\{\mathbf{R}\_{\rm eq}P\_{\rm r0} + \mathbf{X}\_{\rm eq}Q\_{\rm r0}\}}\tag{4.6}$$

$$\frac{\partial V\_{\rm r}}{\partial Q\_{\rm r}} \approx -\frac{\left(\mathbf{X}\_{\rm eq}^{2} + \mathbf{R}\_{\rm eq}^{2}\right)\mathbf{P}\_{\rm r0} + \mathbf{X}\_{\rm eq}}{1 + 2\left(\mathbf{R}\_{\rm eq}\mathbf{P}\_{\rm r0} + \mathbf{X}\_{\rm eq}Q\_{\rm r0}\right)}.\tag{4.7}$$

It can be seen from Eq. (4.6) and (4.7), that in case of <sup>X</sup>eq Req ≫ 1, as in high voltage transmission grids, the voltage sensitivity to reactive power is much higher than the voltage sensitivity to active power, and therefore, the voltage variations depends mainly on the changes in the reactive power. However, for a low <sup>X</sup>eq Req ratio, as common in low voltage distribution grids, the active power variations have also a significant impact on the voltage, since the voltage sensitivity to active power is also high. In fact, for the case of <sup>X</sup>eq Req ≈ 1, the active power has as much influence on the voltage as does the reactive power, since Eq. (4.6) and (4.7) become identical. The difference in voltage sensitivity at different voltage levels is shown more intuitively in Figure 4.2, where the line gradients in Figure 4.2(a) and (b) determine the voltage sensitivity coefficients in high voltage transmission grids and low voltage distribution grids, respectively.

Figure 4.2: Different voltage sensitivity coefficients in (a) high voltage transmission grids, and (b) low voltage distribution grids (figures adopted from [170]).

For a large interconnected power system with numerous buses and lines, the voltage sensitivity coefficients are derived using one of the main three classical approaches described in [171], which all use the load flow calculations. The most common method, which uses the Newton-Raphson load flow formulation is described briefly in the following.

Assuming a balanced grid, by linearizing the load flow equations around one operating point we have

$$
\begin{bmatrix}
\underbrace{J\_{P\mathcal{S}}}\_{J} & \frac{J\_{PV}}{J\_{QV}}
\end{bmatrix}
\begin{bmatrix}
\Delta \mathcal{S} \\
\Delta V
\end{bmatrix} = \begin{bmatrix}
\Delta P \\
\Delta Q
\end{bmatrix}.\tag{4.8}
$$

in which Δ and Δ are small variations in the active and reactive power around the operating point, Δ and Δ are the resulted changes in the voltage phase and magnitude, respectively, and is the Jacobian matrix of the grid. By multiplying both sides of Eq. (4.8) by the inverse of the Jacobian matrix, we have

$$
\begin{bmatrix}
\Delta\delta\\\Delta V
\end{bmatrix} = \underbrace{\begin{bmatrix}
J\_{P\delta} & J\_{PV} \\
\hline
J\_{Q\delta} & J\_{QV}
\end{bmatrix}}\_{\mathcal{S}}^{-1} \begin{bmatrix}
\Delta P\\\Delta Q
\end{bmatrix}.\tag{4.9}
$$

We define the inverse of the Jacobian matrix, in Eq. (4.9), as the voltage sensitivity matrix of the grid. The elements of the voltage sensitivity matrix can be grouped into four submatrices as shown in Eq. (4.10).

$$S = \begin{bmatrix} S\_{\delta P} & S\_{VP} \\ S\_{\delta Q} & S\_{VQ} \end{bmatrix} = \begin{bmatrix} \frac{\partial \delta}{\partial P} & \frac{\partial V}{\partial P} \\ \frac{\partial \delta}{\partial \delta} & \frac{\partial V}{\partial V} \end{bmatrix} \approx \begin{bmatrix} \frac{\Delta \delta}{\Delta P} & \frac{\Delta V}{\Delta P} \\ \frac{\Delta P}{\Delta Q} & \frac{\Delta V}{\Delta Q} \end{bmatrix} \tag{4.10}$$

In Eq. (4.10), corresponds to the submatrix of , which is the partial derivative of the voltage amplitude at bus with respect to the active power changes at the same bus, i.e., the voltage sensitivity to active power changes at the bus , and

$$S\_{VP} = \begin{bmatrix} \frac{\partial V\_2}{\partial P\_2} & \cdots & \frac{\partial V\_2}{\partial P\_n} \\ \vdots & \ddots & \vdots \\ \frac{\partial V\_n}{\partial P\_2} & \cdots & \frac{\partial V\_n}{\partial P\_n} \end{bmatrix} \tag{4.11}$$

Similarly, corresponds to the elements of or the voltage sensitivity to the reactive power variations. The voltage phase angles and their sensitivities are generally not a concern at the low voltage level.

Up to this point, a balanced operation of the grid was assumed for the voltage sensitivity calculations. However, low voltage distribution grids often have asymmetric currents and voltages caused by single-phase loads and generation units. For an unbalanced threephase system, the voltage sensitivity coefficients must be calculated for each phase individually with respect to the power flow of all three phases. Therefore, for bus ,

(4.12)


A more detailed analysis of voltage sensitivity calculations in unbalanced grids can be found in [171] and [172].

## **4.1.2 Limitations of the Classical Method for Calculating the Voltage Sensitivity**

As shown, the classical methods for estimating the voltage sensitivity coefficients use the linearized load flow equations around a single operating point. This approach can provide accurate results only at the particular operating point, in which the calculations are conducted, and neglects the system nonlinearities. However, the grid's operating point varies significantly with time. It has been shown recently that the variations in the voltage sensitivity coefficients can be significant with the changes in the system operating point during a single day [172], [173]. Therefore, accurate estimation of the voltage sensitivity coefficients in every operating point is necessary, in particular when these values are used actively in controller actions, as in [169].

But most importantly, the classical method requires a validated and accurate grid model, or at least its Jacobian matrix, which is not always available. Grid models at the low voltage level are less common, in comparison to higher voltage levels. This can be due to the larger number of these grids in a single region, and the fact that they were considered traditionally the passive parts of the grid. Even if such models are available, they are often not validated nor continuously updated, and are often designed for the planning studies conducted by the grid operator for the expected load in the future. In addition, inverter-based distributed energy resources such as PV systems and their controllers have an impact on the aggregated load/voltage sensitivity [110], but are often not represented accurately in these load flow models.

Due to the aforementioned limitations of the classical approaches, data-driven methods for estimating the voltage sensitivity coefficients and without any knowledge of the grid model can be very advantageous, in particular for low voltage grids. There have been only two attempts for such solutions previously published in literature. Authors in [174], propose an algorithm for estimating the voltage sensitivity coefficients based on the smart meter data. Nonetheless, this method requires simultaneous measurements at different locations in the grid, and historical measurement data and load profiles to train the proposed algorithm. Another work proposes a near real-time estimation method to calculate these coefficients [172]. This work also uses a linear approximation during the measurement period for the estimation procedure. It should be noted that neither of the proposed approaches has been used for allocation studies.

In this thesis, a new data-driven approach has been proposed based on the concept of mutual information. Mutual information has the capability to capture and quantify the nonlinear dependencies between two variables, and it can also generate an output in real time [175]. It is suggested in this thesis that the mutual information can reflect the relative voltage sensitivity coefficients among different locations in a grid, and it only requires measurements at the points of interest. This can be used for applications such as allocation of distributed energy resources such as the allocations of the FESS for power smoothing applications, as demonstrated later in this chapter. Mutual information is described in the following subsection.

#### **4.1.3 Introduction to Mutual Information (MI)**

Mutual Information (MI), originated from the information theory [26], can capture and quantify the dependencies between two variables. It attempts to quantify "Given the values of one variable, how much does it help to guess the values of the other variable?". It reflects the expected amount of information that one variable can reveal about the value of another variable. In our use case, this means how much information about the variations of the active power drawn from a certain point in the grid, can be exposed by just looking at the voltage measurements. Mutual information is much more general than the Pearson correlation coefficient, since it can capture both linear and nonlinear dependencies [176]. Mutual information is zero, if and only if the variables are statistically independent. This is not true for the correlation coefficients.

Because of this ability to capture such complex dependencies, mutual information has been applied in different disciplines, such as machine learning [177] and medical imaging [178]. Mutual information has a great potential to be used in power systems, as highly nonlinear systems. However, to our knowledge, no application of mutual information in power systems has been reported, and this work is the first one to propose it.

Since mutual information is defined originally for probabilities, it requires some estimation to be used on measurement data, which are always limited in the number of data points. A well-known estimate, referred to as the KSG estimate after its creators, Kraskov, Stögbauer, and Grassberger [179], overcomes this challenge by using the distances between similar measurements, rather than estimating the distribution densities. The KSG estimate is fast to compute and yields good results even for small and noisy measurement data [180].

Figure 4.3 presents an example of measurement points for defining the notations required for the KSG estimate of MI. Let = {m<sup>1</sup> , … , m<sup>n</sup> }, in which m = (a , a ) be a collection of measurements of two variables, which are the voltage and the power in phase a in this use case. In Figure 4.3, consists of 13 measurement points. For any ∈ ℕ <sup>+</sup>, the *th* nearest neighbor of a measurement point m is represented by (m ), using the matrix distance of (m , m) = max(|a −a |, |a − a |). For instance, in Figure 4.3, the measurement point m<sup>3</sup> is the fourth nearest neighbor of m<sup>5</sup> , in other words, m<sup>3</sup> = 4(m<sup>5</sup> ). Furthermore, <sup>a</sup> (m ) and <sup>a</sup> (m ), referred to as the marginal counts, are the number of measurement points of voltage and power, respectively, that are not further from m than the (m ). For example, in Figure 4.3, <sup>a</sup> (m<sup>5</sup> ) = 3, i.e., there are three measurement point of <sup>a</sup> that their distance is not larger than the distance of 4(m<sup>5</sup> ) or m<sup>3</sup> . In a similar way, it can be seen that <sup>a</sup> (m<sup>5</sup> ) = 7.

Figure 4.3: Examples of notations used for estimating the mutual information using the KSG estimate.

With these definitions, the mutual information using the KSG estimate for these set of measurements is calculated as [179]

$$MI(V\_{\mathfrak{a}}, P\_{\mathfrak{a}}) = \lambda(n) + \lambda(k) - \frac{1}{k} - \frac{1}{k} \sum\_{l=1}^{n} \left( \lambda\left(c\_{V\_{\mathfrak{a}}}(m\_l)\right) + \lambda\left(c\_{P\_{\mathfrak{a}}}(m\_l)\right) \right). \tag{4.13}$$

in which,

$$
\lambda(\mathbf{x}) = \left(\sum\_{m=1}^{\mathbf{x}-1} \frac{1}{m}\right) - \mathbf{y}.\tag{4.14}
$$

In Eq. (4.13) and (4.14), is chosen to be 4, as recommended in literature [180], and equals to -0.577, which is the Euler-Mascheroni constant.

The value calculated using Eq. (4.13) is referred to as the mutual information score. It is proposed in this work that the mutual information score between voltage and active power measurements can reflect the relative voltage sensitivity to active power changes at different locations in a grid. Similarly, the mutual information score between the voltage and the reactive power can show the relative voltage sensitivity to reactive power variations. This is demonstrated in the next subsections, using real measurement data and the use case of allocation of a FESS in a low voltage distribution grid.

### **4.1.4 The Proposed Allocation Method for a FESS**

In this section, a real case of allocation of a FESS for power smoothing applications in a grid in the southern Germany is presented, as a use case of the proposed method for estimating the relative voltage sensitivity using the mutual information scores.

The distribution system operator intends to install a high-speed FESS in its grid for improving the power quality. The system operator has limited the installation locations to the low voltage side of all the 10/0.4 kV substations, where there is sufficient utilityowned land and electrical and communication infrastructure for the FESS. This results in a total search space of over 1200 substations, as possible candidates for the installation of the FESS. The main application of the FESS is considered to be smoothing out the power drawn from the medium voltage grid in order to reduce the voltage deviations. Therefore, a location with a high voltage sensitivity to active power changes is desired.

To find the substation with the highest voltage sensitivity, a new allocation algorithm based on the concept of mutual information for estimating the relative voltage sensitivity is proposed. The proposed allocation algorithm consists of two main steps, which are described in detail in the following, along with the data collection procedure. The flow chart of the proposed allocation algorithm is shown in Figure 4.4.

#### **4.1.4.1 Step 1: Selecting the Top Candidates for Measurements**

Medium voltage grids are often monitored only at their interconnection with the high voltage grids, while low voltage grids are usually not monitored at all [181]. This is commonly justified by high expenses or being unnecessary, since they are large infrastructures, and were assumed to be the passive parts of the grid. Many distribution grid operators, collect measurement data only, if power quality issues or interruptions are reported. Despite the ongoing effort to improve the observability at these voltage levels in order to deal with the increasing share of distributed energy resources, the number of fixed measurement devices at these voltage levels is still low. This has been also the case for the medium voltage and low voltage grids investigated in this work. Therefore, first, measurement data have to be collected.

Figure 4.4: The flow chart of the proposed allocation algorithm for the FESS based on mutual information.

Since collecting measurement data from all 10/0.4 kV substations is not feasible nor practical, the substations with the highest potential for the FESS installation are selected for conducting measurements. As measurement data are not available at this point, the voltage sensitivity coefficients are calculated using the classical method based on the Newton-Raphson load flow calculations at the nominal operating point, as described in subsection 4.1.1. A static load flow model of the medium voltage grid in DIgSILENT PowerFactory is provided by the system operator, which is used for these calculations. It should be noted that this model is used only for planning studies of the grid operator.

This first step of the allocation algorithm is programmed using the DIgSILENT Programming Language (DPL) in DIgSILENT PowerFactory. Firstly, all the 10/0.4 kV substations are selected and the voltage sensitivity coefficients are calculated at the nominal operation point at all buses, using the method explained in subsection 4.1.1. The top nine buses with the highest voltage sensitivity to active power according to the model

have been chosen as the most attractive candidates for conducting measurements and possibly the installation of the FESS. Next, the FESS is installed at each of these buses in the model, and the load flow calculations are repeated to check the following technical constraints (C1-C4):


In order to carry out the load flow analysis after the FESS installation, the bus at which the FESS is connected, is modelled as a PQ bus, in which the reactive power follows the Q(U) characteristics given in the latest German grid code for the low voltage grids, i.e., the VDE-AR-N 4105:2018-11 [182]. The Q(U) characteristic for the ESS connected to the low voltage grids is illustrated in Figure 4.5(b), together with the P(f) characteristic. In this figure, <sup>g</sup> and <sup>g</sup> are the voltage and frequency of the grid, respectively, max is the maximum reactive power of the ESS, and Emax is its maximum active power. Different load flow models for an ESS can be found in [183].

Figure 4.5: Active and reactive power control characteristics according to the German grid code VDE-AR-N 4105:2018-11. (a) The P(f) characteristic for frequency deviations, and (b) the Q(U) characteristic for voltage deviations.

If one or more of these constraints are violated after the installation of the FESS at a particular bus, the bus is replaced with the next bus with the highest voltage sensitivity to active power changes, and the constraints are checked again. There are also some nontechnical constraints, such as having adequate space and infrastructure for installing the FESS. If the selected location does not have any of these conditions, the location is also replaced with the next most sensitive bus. The top nine candidates which fulfil both the technical and non-technical requirements are chosen for field measurements.

#### **4.1.4.2 Measurement Setup at the Low Voltage Distribution Networks**

In this subsection, the measurement setup used at the most attractive candidates for the FESS installation in a grid in southern Germany are described.

The measurements were conducted in a close cooperation with the distribution system operator of the region. An exemplary schematic of the measurement setup is illustrated in Figure 4.6. For each measurement, the measurement device was installed at the secondary side of the 10/0.4 kV transformer, in order to measure the power profile of the whole low voltage distribution grid, including all its feeders. Voltages were measured directly, while currents were measured using 800/5 A current transducers. The neutral current was measured separately.

Figure 4.6: Schematic example of a measurement setup at the low voltage distribution networks.

The measurements were conducted using A. Eberle's power quality and disturbance recorder, "PQI-DA smart" [184]. This device samples the data at the high resolution of 41 kHz, but records the data on longer intervals, starting from every second. The highresolution data is only recorded, if one of the predefined thresholds, determined by the power quality indices of the European standard DIN EN 50160 [159], are violated. When a continuous recording at a high resolution was required, the in-house built Electrical Data Recorder (EDR) [185] was used, which is capable of measuring and recording continuously at 25 kHz of sampling rate.

The measurement data were collected at the nine substations, which were chosen according to the step one of the allocation algorithm. A total of 223 days and 124 GB of data has been recorded. The measurements were conducted on low-voltage grids with different types of loads, including residential, commercial, industrial, and a combination of these types. Therefore, different patterns of load profiles were observed. Also, the low voltage grids have shown different penetration levels of PV generation, from no PV installations at all, to the case with enough PV generation to surpass the local demand at peak of generation. A minimum of two weeks of measurements have been conducted at the most attractive locations for the FESS installation.

From the perspective of power quality, no major violations of the European standard DIN EN 50160 [159] has been observed during the whole measurements. However, a series of interesting events have been recorded, which are presented in Appendix A.

#### **4.1.4.3 Step 2: Estimating the Mutual Information Scores**

For each bus, the mutual information scores between the voltage and active and reactive power measurements are calculated according to Eq. (4.13) using a window length of two weeks. The average of the mutual information scores of the three phases at the top nine candidates is shown in Figure 4.7.

As seen in Figure 4.7, bus 4, followed by bus 3 and 1, shows the highest mutual information scores between the voltage and the active and reactive powers, indicating a higher dependency of the voltage and the power at these points of the grid. This shows a higher voltage sensitivity to active and reactive power changes, in comparison to the other buses. Therefore, they can be more suitable for the installation of the FESS.

Bus 5 and 2 have the lowest mutual information scores, indicating a low voltage sensitivity to the active and reactive power changes. This can imply a strong grid connection, where active and reactive power variations have less impact on the voltage, when compared to other candidates. In addition, the results for bus 5, 7, and 8 show a higher mutual information score between voltage and active power, in comparison to the voltage and reactive power. This reflects a higher voltage sensitivity to active power in comparison to the reactive power, which is the characteristic of a low X/R ratio.

Figure 4.7: Comparison of the mutual information (MI) scores of the top nine candidates for the FESS installation.

#### **4.1.5 Evaluation of the Proposed Allocation Method**

In order to evaluate that the mutual information scores can reflect the relative voltage sensitivity at the different locations in a grid, the voltage and active power measurements, collected at bus 4 and 5, i.e., the buses with the highest and lowest mutual information scores between voltage and active power, are illustrated in Figure 4.8. Measurement data are shown for a period of two days with a 1-minute data resolution. The measurements are shown for bus 4 and 5 specifically, as their differences can be seen more clearly, compared to two buses with similar mutual information scores.

By observing Figures 4.8(a) and (b), it can be seen that the active power changes have a clear and significant impact on the voltage at bus 4, and these two variables greatly depend on each other. It can be observed, that the voltage variations follow the changes in the active power in the opposite direction, and their general pattern seems very similar. This indicates a weak connection point of the grid. For bus 5 with the lowest mutual information score, it can be observed from Figures 4.8(c) and (d), that the effect of the active power variations on the voltage is not easily recognizable. The dependency between the voltage and active power seems to be much lower, when compared to bus 4, and the patterns in the active power changes cannot be easily observed in the voltage. This is the characteristic of a stronger connection point of the grid, in comparison to bus 4. Similar results can be observed by comparing other buses to each other and by comparing the sensitivity coefficients of reactive power. Therefore, it can be concluded the proposed data-driven method based on the concept of mutual information can reflect

the relative voltage sensitivity to active and reactive power changes at different locations in the grid.

Figure 4.8: Measurement data at bus 4 and 5 illustrated for a period of two days with a 1-minute data resolution. (a) Bus 4 voltage, (b) bus 4 active power, (c) bus 5 voltage, and (d) bus 5 active power. Bus 4 has the highest mutual information score between voltage and active power, while bus 5 has the lowest score.

The proposed approach can be used for a variety of applications such as the allocation of distributed energy resources and ESS, in particular for improving the voltage quality. For the presented use case of allocation of a FESS in a grid in southern Germany, the proposed allocation algorithm has shown that the bus 4, followed by bus 3 and 1, has the highest potential for installing the FESS, as a higher relative voltage sensitivity to active power variations is observed (see Figure 4.8). However, just having a high voltage sensitivity to active power changes is not adequate in order to effectively benefit from a FESS. It is also important to have short-term active power fluctuations with high peaks, such as the ones caused by pulsed power loads, at the selected location in the grid. Such peaks in locations with a high voltage sensitivity to active changes can lead to significant fast voltage changes, which is a power quality violation. Where there is a high voltage sensitivity to active power variations, but no sharp active power fluctuations are present, significant fast voltage changes caused by active power changes will not be observed.

Therefore, the measurement results at the top three candidates with the highest mutual information scores, bus 4, 3 and 1, have been compared in terms of having short-term active power variations. Measurements result show, that the active power at bus 1 shows the highest number of short-term active power peaks. Due to the high voltage sensitivity at this location, the highest number of voltage deviations above 3 % is also observed at this location during the whole measurement period. A total number of 76 incidents of voltage deviations above ±3% of the nominal value were recorded during the measurement period of two weeks, which were caused by sharp active power changes. There are fewer voltage deviations recorded at bus 3 and 4, and none of them were caused by active power changes. An example of variations of the voltages and currents during several hours of measurements at bus 1 is shown in Figure 4.9. It is clear from Figure 4.9 that a pulsed power load is connected at this point in the grid. Since the voltage sensitivity to active power changes is also relatively high at this bus, these pulsed power loads cause significant fast voltage changes, as seen in Figure 4.9.

In conclusion, bus 1 is chosen as the most attractive location for the FESS installation among more than 1200 substations in the studied German grid. The grid can benefit the most from this storage technology at this particular location, as it is a weak point in the grid, and a large pulsed power load is also connected to it. The FESS can smoothen the power at this substation, by providing the short-term power requirements. Subsequently, the voltage deviations caused by the power pulsed load can be avoided.

Figure 4.9: An example of the variations in the voltages and currents at bus 1, the most attractive candidate for the FESS installation during a period of 6 hours.

# **4.2 A Data-driven Method for Sizing ESS using Standard Patterns<sup>2</sup>**

A high penetration level of PV generation in low voltage distribution grids can cause several challenges for the distribution system operator. The peak in PV generation appears around noon, when the demand is often low in residential areas with a high fulltime employment rate [186]. This can lead to a reverse power flow from the PV installations towards the medium voltage grid during noon, which can potentially lead to over-voltage issues in the low voltage grids, due to the low X/R ratios [187], [188]. This can limit the PV hosting capacity of the low voltage grid for new PV installations [169], [187]. An example of a daily power profile of a low voltage distribution grid with a high share of PV generation is shown in Figure 4.10. This data is collected at a 10/0.4 kV substation in southern Germany during a summer day in 2018, using the setup described in subsection 4.1.4.2. It can be seen, that the active power drawn from the substation becomes negative after around 11:00 until approximately 15:00. The problem of reverse power flow and non-coinciding demand and generation can be solved by installing an ESS with high energy density, such as a lithium-ion battery energy storage system, if sized properly.

Figure 4.10: Example of a daily power profile of a low voltage distribution grid with a high share of PV generation during a summer day of 2018 in southern Germany.

In addition, as discussed in the previous section, the variations in the active power can significantly impact the voltage in grids with a low X/R ratio. Therefore, the fast variations in the solar irradiance, e.g., caused by passing clouds, which cause the power

<sup>2</sup> This section presents an extended verion of the following publication: S. Karrari, N. Ludwig, M. Noe, V. Hagenmeyer, "Sizing Centralized Energy Storage Systems in Distribution Networks Using Standard Patterns," in IEEE PowerTech 2019, Milan, 2019.

output of the PV systems to fluctuate, or variable industrial loads, can cause rapid voltage deviations in such grids. These rapid voltage fluctuations might violate the power quality threshold for the maximum permissible fast voltage changes [189], such as the 0.05 p.u. limit, given in the European grid code DIN EN 50160 [159]. Examples of large active power fluctuations can be observed in the measured power profile of Figure 4.10, between 13:00 and 15:00. This fast active power changes can be compensated using a high power density ESS with a large number of cycles such as a FESS, if it has the sufficient rated power and capacity.

The objective of this section is to derive the energy storage characteristics from the recorded data. As seen in the related work, the first step for sizing an ESS based on historical measurement data is to select a suitable data set as the input for the sizing methodology. The data resolution and the chosen period of the recorded data can both impact the sizing outcome. While 15-minute [190] or even hourly data [191] are frequently used in storage sizing studies, it is recommended to use higher resolution data, such as a data resolution of 1 second, in order to study the use of an ESS in the presence of intermittent PV generation in low voltage grids [158]. Similarly, another study concluded that the use of the full frequency spectrum for storage sizing studies is highly recommended [192]. However, using high-resolution data, or even low-resolution data on a long horizon (several years) leads to large data sets. This can drastically increase the computation time required for solving the sizing problem, in particular for complex nonlinear optimization frameworks, which can make the problem intractable [193]–[195]. In order to avoid dealing with such large data sets, it is common practice to choose some specific short period, e.g., one day, as a representative of the measurement data for the sizing study, or to reduce the data resolution by aggregating the data. However, the choice of a representative day is often not clearly explained, nor properly justified. Frequently, an arbitrary day or some worst-case scenario, e.g., a day with high PV generation and low demand, or a day with low PV generation and high load, is selected. These approaches can potentially lead to over- or under-sized ESS. For instance, in [196], a random day is chosen as the load profile for the sizing study, which can lead to very specific storage requirements that are not suitable for other days. Due to the lack of highresolution data, 1-minute data are artificially generated from hourly data in [197], but again, the data of a typical day is used, and it is not clear how this typical day is selected. A similar approach has been reported in [198], where the logic behind choosing the typical day is not provided by the authors. In [192], a random selection of days is used to calculate the storage ratings. The alternative solution, which is reducing the data resolution is also not recommended, as discussed earlier. As an example, authors in [199] formulate the sizing problem as a convex optimization problem, but convert the 10 second load data to 15-minute data to ease computation over a year, which eliminates the high-frequency components of the power measurements.

To properly address this issue in sizing studies, clustering techniques such as K-means clustering [194], [200], [201] and fuzzy C-means clustering [202] were used to group similar daily power profiles and generate representative ones, that can be used for the sizing problem, instead of the whole data set. This is a significant improvement to the arbitrary choice of a typical day. However, when using clustering techniques, it is required to specify beforehand, how many clusters are needed to be found. Depending on the number of clusters and selected features, this can result in power profiles being grouped into one cluster, while not being necessarily similar, since often simple features such as the peak or mean value is used for the clustering (see [202] as an example). These approaches can also generate clusters of power profiles that occur very seldom, which cannot represent the data set properly. It should also be noted that the related works which use clustering techniques all use hourly data for long term scheduling problems and high-resolution data is not used. Table 4.1 provides a summary of the methods used in literature for reducing the data size for sizing studies and their drawbacks.


Table 4.1: Summary of the previously used methods for reducing data size for sizing studies and their drawbacks and limitations.

In this thesis, a novel alternative is proposed, where reoccurring daily consumptions patterns are detected using the motif discovery algorithm [27], which are then used for the sizing of the ESS, rather than using an arbitrary day, clustering techniques or reducing the data resolution. Motif discovery uses the Symbolic Aggregated Approximation (SAX) of the data set, instead of the raw data, which helps to find similar, but yet, not the same patterns that are being repeated during the whole data set. It also uses the random projection algorithm [204], which makes the process of finding similar patterns faster than the clustering approaches [205]. Moreover, motif discovery only groups power profiles when they are very similar. There is also no need to specify the number of patterns that needs to be found, and patterns or *motifs* are generated only if they are being repeated often. The pattern with a high number of reoccurrences and high probability are then used for deriving the characteristics of the ESS, rather than using some random day or the worst-case scenario. Since the length of pattern used for the sizing data is only one day, high-resolution data can be easily employed with no concern on the computational burden of the sizing algorithm.

In this work, the motif discovery algorithm has been applied for sizing two different types of ESS at 10/0.4 kV substations for solving different challenges in the low voltage distribution grids:


For deriving the characteristics of each ESS, firstly, reoccurring daily consumption patterns are detected using the motif discovery algorithm from the measurement data collected at several low voltage grids in southern Germany. Then, a low-pass filter is applied to decouple the discovered patterns into high- and low-frequency components for sizing different storage technologies. The cut-off frequency of the low-pass filter is determined through a simple optimization framework. The nominal capacity, nominal power, maximum ramp rate, and the number of times per day where each ESS changes from charging to discharging mode are calculated based on the power profiles allocated to each type of storage technology. The proposed sizing methodology is evaluated by assessing the performance of both ESS types, with the characteristics derived from the standard patterns, during the time horizon of the whole data set. In addition, the effects of data resolution, the cut-off frequency of the low-pass filter, and choosing patterns other than the most reoccurring one on the sizing outcome are also investigated. Although the proposed method is used for sizing a centralized ESS at the 10/0.4 kV substations, this method can be applied for sizing all types of ESS for any other application using historical measurement data. Examples include behind-the-meter applications of an ESS, such as increasing the self-consumption of the PV generated power or energy arbitrage.

The rest of this section is organized as follows. The collected measurement data used for deriving the storage characteristics in this work are presented in the next subsection, followed by a description of the motif discovery algorithm. The discovered standard patterns and their characteristics are presented in the following subsection. Next, the most important storage characteristics are formulated using the power profile that each ESS has to cover. The procedure to determine the optimum cut-off frequency for applying a low-pass filter on each pattern is described next with the aim optimizing the characteristics of the two major types of ESS. The sizing outcome and the effects of data resolution, cut-off frequency of the filter, and selecting other daily patterns on the storage characteristics are discussed in the following subsection. Lastly, the evaluation of the proposed sizing methodology is presented using the original full data set.

#### **4.2.1 The Input Data for the Sizing**

This work aims for installing a centralized ESS at a 10/0.4 kV substation. Therefore, the input data for the sizing study are collected at several 10/0.4 kV substations in southern Germany in the summer and fall of 2018, using the setup described in subsection 4.1.4.2. In comparison to the distributed ESS, a centralized ESS at the substation has the advantage of a single grid interface, a simpler control design, access to the electrical and monitoring systems of the substation, and availability of utility-owned land [206].

Among the substations from which data have been collected, four substations were selected for the sizing study, which represent different levels of PV penetration, and different types of loads. Moreover, the measured data has been recorded using three different resolutions of 1 second, 1 minute, and 10 minutes. Therefore, the effects of data resolution on the sizing of different storage technologies can also be investigated. A minimum of two weeks of measurement data was recorded at each substation.

As an example, the collected data with 1-minute resolution during a period of 8 days, starting from a Monday, has been illustrated in Figure 4.11 for the four selected buses. Clearly, there are patterns that are repeating themselves on a daily basis at each bus. While these patterns seem quite similar, they are also not exactly the same. For bus 1 and 2, which their power profiles are shown in Figures 4.11(a) and (b), there are different patterns occurring during the weekends (the 6th and 7th day). This is due the fact that these low voltage grids have mostly industrial and commercial loads. However, starting from the day 8, the pattern for a working day starts to reappear. In bus 1, the share of PV generation in comparison to the power consumption is so low that it can hardly be observed in the power profile during the working days, and it can only be seen during the weekends. In the power profile of bus 2, a slightly higher share of PV generation is observed, in which the increase in the consumption around noon during working days is more or less covered by the PV generation at its peak. In both cases of bus 1 and 2, negative power flow is observed only at weekends, when the load is at its minimum. Bus 3 and 4, which their power profiles are shown in Figures 4.11(c) and (d), represent two residential low voltage distribution grids with a relatively high level of PV penetration. In these cases, the consumption patterns are relatively similar for both workdays and weekends. The power profile measured at bus 3 shows a high peak only during the evening, while in the power profile of bus 4, the morning peak is as large as the evening peak. Moving from bus 1 towards bus 4, the ratio of PV generation to the peak consumption power increases, up to a point that in case of bus 3 and 4, the PV generation is dominant and leads to reverse power flows during noon in several days. Large power fluctuations around noon are also more present in buses with a higher share of PV generation, such as in bus 3, which could have been caused by passing clouds over the PV panels. It is important to note that these measurements were not conducted simultaneously, as a single measurement device was used.

Figure 4.11: Active power measurements at (a) bus 1, (b) bus 2, (c) bus 3, and (d) bus 4. The data is shown with the resolution of 1 minute for a period of 8 days, starting from Monday, which makes the 6 th and 7th days, weekends.

The power profiles at different buses also differ in terms of their frequency components. This can be illustrated using the Discrete Fourier Transform (DFT) applied to the 1 seccond measurement data. The one-sided frequency spectrum using the DFT is shown in Figure 4.12 for each bus. Each frequency component represents in this figure a periodic cycling of power. As seen, each power profile has a different amplitude at each frequency. For instance, the frequency component of 0.0116 mHz, which corresponds to every 24 hours, can be an indication of daily reoccurrence of a consumption pattern, which is highest for bus 4, followed by bus 3. This is due to having similar power profiles during the whole week, including weekends. The power profile at bus 4 shows a higher magnitude at 0.0231 mHz in comparison to other buses, which corresponds to 12 hours. This can clearly be seen from Figure 4.11(d) as well, where each day can be divided into two similar patterns. High frequency components above 2 mHz are observed mostly at bus 1, as it consists of mostly industrial loads with a high power variability.

Figure 4.12: One-sided frequency spectrum of the power profiles at all buses using the 1-second data for a period of one week.

#### **4.2.2 Introduction to Motif Discovery**

As previously discussed, the first step towards sizing an ESS using historical measurement data is to select the power profile, which can accurately represent the whole data set. Instead of using an arbitrary day or the worst-case scenario, which can lead to very specific sizing outcomes, we detect and use a reoccurring daily consumption pattern at each substation using the motif discovery algorithm. The approximation of the whole time series with a standard reoccurring pattern helps us to focus on regular reoccurring patterns in the consumption, rather than on single occurrences.

Finding reoccurring patterns in time series is a common unsupervised learning problem. In this work, we propose the use of motif discovery, introduced in [27], and its improved version in [207], tailored for energy-related time series. A major advantage of using motif discovery is the use of the Symbolic Aggregated Approximation (SAX) of the time series rather than the time series itself for finding similar patterns. This allows to find patterns that are similar and yet, not the same. In contrast to the established clustering approaches [194], [200]–[202], there is also no need to specify beforehand the number of groups or clusters, and power profiles are only grouped together, when they are found to be similar. Power profiles that are rarely reappearing are automatically left outside.

Motif discovery has been applied in various disciplines for a variety of applications, including identifying seismic activity patterns in active volcanos [208] and classifying heartbeat sounds [209]. In the realm of energy, motif discovery has been applied for finding flexibility potentials in industrial processes [210] and non-intrusive load monitoring [211]. Prior to this work, there are no reports on using motif discovery for storage sizing studies or other applications in power systems.

The algorithm consists of five main steps, as shown in Figure 4.13, which are described in the following according to [207]. The term *motif* refers to a pattern or a series of sequences in the time series that are repeated often.

Figure 4.13: The process of the motif discovery using SAX and random projection (adopted from [207]).

**Sequencing algorithms:** In the first step, the whole time series data, which consists of data points, is divided into sequences with equal length of . As we are looking for daily patterns in this work, the length of the sequence is set to correspond to 24 hours, depending on the data resolution. For instance, for the 1-second data, is 86400.

**Symbolic Aggregated Approximation (SAX):** Next, the sequences are converted into a series of letters out of the alphabet, similar to the DNA sequences, from which this method was inspired. To do so, the sequence time series Y is divided into equal pieces with a predefined length of . All data points within each piece are aggregated by their mean value, <sup>i</sup> PAA. This is known as Piecewise Aggregated Approximation (PAA) of the data, which helps to reduce the effect of noise and minor variations in finding patterns.

In the following step, a letter out of the alphabet is assigned to each piece, according the value of <sup>i</sup> PAA. First, a section of the alphabet with the predefined length of together with a list of breaking points are defined, in which = <sup>1</sup> , <sup>2</sup> , … , a−1 , and i−1 < <sup>i</sup> . The breaking points are chosen to be the quantiles of the empirical cumulative distribution function of the sequence, which divides the area underneath into equal areas. We assign = 10, so that each letter corresponds to a meaningful quantile of the distribution function. If the mean value of each peace (<sup>i</sup> PAA) is in the range of [j−1 , j), then its SAX representation (<sup>i</sup> SAX) is the corresponding letter of <sup>j</sup> , or

$$\mathbf{y}\_{\parallel}^{\text{SAX}} = \alpha\_{\parallel} \qquad \text{if, } \mathbf{y}\_{\parallel}^{\text{PA}} \in \left[\mathcal{J}\_{\parallel -1}, \mathcal{J}\_{\parallel}\right]. \tag{4.15}$$

For observations = (<sup>1</sup> , <sup>2</sup> , … , ), the empirical cumulative distribution function () is the fraction of observations which are less than or equal to , i.e.,

$$F\_n(p) = \frac{1}{n} \sum\_{l=1}^n I(\mathbf{x}\_l \le p). \tag{4.16}$$

At the end of this step, we have strings of letters or *words* representing each sequence, which in our use case are the daily power profiles. The PAA and SAX of a time series are illustrated more intuitively in Figure 4.14. In this figure, the original time series is shown in grey, the black line shows the PAA of the data, and the blue letters show its SAX representation.

Figure 4.14: An example of PAA of a time series (black line) and its SAX representation (the blue letters). The time series itself is shown with the grey line (adopted from [207]).

**Random projection:** To find the potential motifs, one can simply compare each sequence to all other sequences to look for similarities, but this can be time consuming for large data sets. Instead, motif discovery uses random projection [204], which converges much faster. The SAX representations of the data are saved in a similarity matrix ∈ ( − + 1) × as rows, using a sliding window. Here, is the number of observations, and is the number of sequences. In each iteration of the random projection algorithm, columns of the similarity matrix are selected randomly, where ≤ , and is a user-defined parameter. The words that are built with the selected letters are compared with all rows of . When the words are similar, the corresponding entry in the so-called collision matrix is incremented. The entries with the highest values in the collision matrix are considered potential motifs. The motifs are then transformed back to their original data form, which are the daily power profiles in this use case.

**Candidate evaluation:** The aim of this final step to find how often the candidate motifs occur during the whole time series. The potential motifs are iterated over the original time series to find the instances where the motifs occur, which are the days that follow the same pattern. Therefore, while finding motifs is based on their SAX representation, their evaluation is carried out using the original time series.

Since the length of the sequences are equal in this use case, a simple Euclidean distance measure is used to find the occurrences of each possible motifs. For two sequences of <sup>i</sup> and <sup>j</sup> , (<sup>i</sup> , ) is the Euclidean distance between them, and is calculated according to

$$d\left(y\_{l}, y\_{l}\right) = \sqrt{\Sigma\_{t=1}^{n} \left(y\_{l}(t) - y\_{l}(t)\right)^{2}}.\tag{4.17}$$

Each sequence similar to a motif is considered an occurrence of the motif. Each motif can have several occurrences within the whole time series. We define the *standard pattern* as the 80 %-quantile of all the occurrences of the motif with the highest number of occurrences. This quantile is chosen, as the loads are also supported by the grid in this use case. Choosing a higher percentage quantile is also possible, which will lead to larger storage units.

#### **4.2.3 Standard Patterns for Sizing ESS**

Using the described motif discovery algorithm, the standard consumption patterns at the four low voltage distribution grids are obtained. A minimum of two weeks of measurement data have been used to derive these patterns. However, larger data sets can also be easily used, considering the fast convergence of the random projection algorithm. The standard patterns are derived separately for the 1-second, 1-minute, and 10-minute data for each low voltage distribution grid, in order to study the effect of data aggregation on the sizing outcome.

The standard patterns using the 1-minute data for all buses are shown in Figure 4.15. Similar patterns are derived using the data with other resolutions. As seen in Figure 4.15(a), the standard pattern of bus 1 represents a low voltage grid with smallscale industrial facility, which is active only within working hours, and has a low share of PV generation. The share of PV generation at bus 2 is slightly higher, and just enough to reduce the midday peak to the base load value, while for the buses 3 and 4 the PV generation can almost cover the load for short periods at its peak. The standard pattern at bus 4 shows a slight reverse power flow between 13:00 and 14:00.

Figure 4.15: The standard consumption patterns at the four buses: (a) Bus 1: No PV generation with industrial loads. (b) Bus 2: Limited PV generation to cover only the increase in consumption at peak generation. (c) Bus 3: Significant share of PV generation in a residential grid, resulting in a typical *duck curve*. (d) Bus 4: High share of PV production in a residential grid, leading to slight negative power flow at peak generation.

In order to see the effect of PV generation on the frequency components of the power profiles, a time-frequency analysis on the standard patterns is carried out. Figure 4.16 shows the spectrogram of the standard patterns using the 1-second data of the four buses, which is calculated using the short-time Fourier transform and a Hamming window with the length of 1-minute with no overlaps. The spectrum ends at half of the sampling frequency, since according to the Nyquist principle, the spectrum is symmetrical around this value. As seen in Figures 4.16(c) and (d), at bus 3 and 4 with a high share of PV generation, the frequency content varies significantly with the increase in the PV generation around noon at both high and low frequencies. The constant power lines on the spectrum become vertical from being horizontal, which indicates that the PV generation alters a vast range of frequencies at the same time. With the decrease in PV generation at the evening, the frequency spectrum is restored, and shows more or less the same frequency components until the end of the day.

Another interesting observation is that at bus 1, which supports an industrial facility, both high- and low-frequency components show high values within the activity period of the industrial facility, with the exception of the lunch break at approximately 12:00, as seen in Figure 4.16(a).

Figure 4.16: The spectrogram of the standard patterns using the 1-second data and hamming window of 1 minute length with no overlaps for (a) bus 1, (b) bus 2, (3) bus 3, and (d) bus 4.

The standard patterns of Figure 4.15 are the motifs with the highest number of reoccurrences. However, it is possible that motif discovery algorithm detects several other patterns as well, with lower number of occurrences. This can be the case, for instance, where workdays and weekends follow different consumption patterns, or large data sets such as yearly time series are used, where seasonal power profiles can be quite different. In our data sets, several patterns were detected at bus 2 when using the 1 minute and 10-minute data. When using the 1-minute data, three motifs have been detected at bus 2, which are shown in Figure 4.17. In this figure, pattern 1 represents a normal weekday pattern, which has the highest number of occurrences, and it is the main pattern used for the storage sizing later within this chapter. Patterns 2 shows the consumption pattern during the weekends and public holidays, while patterns 3 represents the power profile of the first working day after a weekend or a holiday. The difference between pattern 1 and 3 can be associated with the loads that stay connected during the workdays but are disconnected for weekends. The effect of using the other patterns with lower number of reoccurrences on the sizing is investigated later in subsection 4.2.6.3, which can also show how an arbitrary choice of days for a sizing study can lead to different sizing outcomes.

Figure 4.17: The three different patterns detected at bus 2 using motif discovery. Pattern 1 is the consumption pattern during a typical workday with the highest number of reoccurrences, while pattern 2 shows the consumption pattern over the weekdays and public holidays. Pattern 3 represents the power profile of the first day after a weekend or a holiday.

#### **4.2.4 Deriving Energy Storage Characteristics**

The standard power profile at each bus is used to calculate the required ESS characteristics. Here, it is assumed that the daily average of the standard patterns is covered by the grid and the ESS provides the variations from the daily average. Using this approach, the evening peak, the reverse power flow, the high ramp rate requirement during the evening hours, and the high-frequency power fluctuations are reduced using the ESS, and the distribution system operator sees a roughly constant demand from the low voltage grid throughout the day.

But different storage technologies operate at different time scales. They also differ substantially in terms of their lifetime and the maximum number of useful cycles. As discussed in Chapter 2, while a FESS can provide a significant number of short-term power cycles with no degradation over its lifetime, a high energy density ESS, such as a lithium-ion battery, are more limited in this regard. Therefore, in our use case, it is suggested to use a high energy density ESS for peak shaving and load balancing, while a high power density ESS such as a FESS is used to cover the short-term power fluctuations. This helps smoothing the charging and discharging profile of the lithium-ion battery, and therefore, preserving its lifetime and capacity [55]. Both flywheels and supercapacitors can be used for smoothing the charging profile of a battery, but it has been recently demonstrated that a hybrid battery-flywheel system outperforms a batterysupercapacitor combination in terms of reducing the peak current drawn from the battery, and the overall system efficiency [57]. Reducing the maximum current drawn from the lithium-ion cells can limit the degradation of the anode active material and lithium plating, as shown in [17].

Therefore, the low-frequency components of the variations in power from its average are allocated to a high energy density ESS with limited number of cycles, referred to as *type 1 ESS*. The high-frequency components are allocated to a high-power density ESS, referred to as *type 2 ESS*, for smoothing out the short-term power variations. Each storage technology can also be used individually, if required. The type 2 ESS, such as a FESS, can be used alone to smoothen the sharp active variations drawn from the grid in order to avoid fast voltage changes in weak parts of the grid, similar to the use case of section 4.1. It can also be used in combination with a type 1 ESS in order to extend its lifetime. The type 1 ESS can also be installed alone to reduce the evening peak, avoid reserved power, and reduce the required ramp rate during evenings. However, in this case, it also has to cover the short-term power fluctuations, which can reduce its effective lifetime. Combining these technologies in a hybrid ESS can provide all these advantages at once.

The following describes how the nominal capacity, the nominal power, the maximum ramp rate, and the number of mode changes from charging to discharging and vice versa, are derived from the allocated power profile to each ESS. These calculations are repeated for each type of ESS, each data resolution, and each standard pattern. Practical aspects such as the energy conversion losses and the requirements to preserve the lifetime of the type 1 ESS are also considered in the sizing calculations.

#### **4.2.4.1 Nominal Capacity**

The nominal capacity of an ESS determines the characteristics of the storage medium, such as the number of cells in a lithium-ion battery, or the inertia and maximum rotational speed in a FESS, as discussed in Chapter 2.

The nominal capacity (En) can be calculated by integrating over the daily power profile allocated to each ESS, and then computing the difference between the highest and lowest peaks of the output signal. The efficiency of the ESS (), which includes the efficiency the power electronics interface and the storage medium, should also be considered, as the ESS should provide higher power than the allocated power while discharging, and lower power is stored in the ESS while charging, considering the conversion losses. The sign function () is used to determine from the power profile that whether the ESS is being charged or discharged for applying the efficiency of the ESS.

In addition, when using storage technologies such as lithium-ion batteries, it can be beneficial to the system's lifetime to avoid very low and high values of state of charge. It has been shown for lithium-ion batteries that avoiding high values of state of charge (above 80%) can significantly help reduce the cathode degradation and calendar aging of the battery [17]. Similarly, deep discharging should also be avoided, as it can significantly increase the internal resistance in the lithium-ion cells. Therefore, it is recommended to reserve a certain non-usable capacity to prevent deep discharging and over-discharging the system [17]. Hence, in this work, the state of charge is kept within the certain range of [SOCmin, SOCmax] in order to preserve the system lifetime. Considering all these factors, the nominal capacity of the ESS is calculated using

$$\mathbf{E}\_{\mathbf{n}} = \frac{ma \, x \Big(\eta^{-\text{sgn}(P(\mathbf{r}))} \int\_0^\mathbf{f} P(\mathbf{r}) d\mathbf{r}\big) - mh \Big(\eta^{-\text{sgn}(P(\mathbf{r}))} \int\_0^\mathbf{f} P(\mathbf{r}) d\mathbf{r}\big)}{\text{SOC}\_{\text{max}} - \text{SOC}\_{\text{min}}}. \tag{4.18}$$

In Eq. (4.18), () is the allocated power profile to the ESS, and corresponds to 24 hours, depending on the data resolution.

For the type 1 ESS, the parameters for a lithium-ion battery are assumed. The minimum state of charge is assumed to be 10 %, while the maximum state of charge is chosen to be 80 %, as recommended in [17]. The roundtrip efficiency is assumed to be 90 %, which is the efficiency of a Tesla PowerWall 2 home battery system [212]. The same efficiency is used for the type 2 ESS, with the minimum and maximum values for the state of charge to be 25 % and 100 %, respectively. This is due to the fact that a FESS have no constraints on being fully charged, but the machine's torque constraint limits the maximum deliverable power at low values of state of charge (see subsection 6.3.1).

#### **4.2.4.2 Nominal Power**

The nominal power (Pn) of an ESS is often limited by its power electronics interface. It is simply the maximum value of the power profile allocated to the ESS, considering both charging and discharging modes, and the efficiency of the ESS. Therefore,

$$P\_{\mathbf{n}} = \max(\left|\eta^{-\text{sgn}(P(\mathbf{t}))}P(\mathbf{t})\right|). \tag{4.19}$$

#### **4.2.4.3 Maximum Ramp Rate**

In this work, the maximum ramp rate (Rn) is considered to be the maximum change in the storage power between two consecutive measurement points in time. Considering the sampling frequency (<sup>s</sup> ) of the measured power signal, the maximum ramp rate is calculated as

$$\mathcal{R}\_{\mathbf{n}} = \max \left( \eta^{-\operatorname{sgn}(p(\mathbf{r}))} P(\mathbf{r}) - \eta^{-\operatorname{sgn}\left(P\left(\mathbf{r} - \frac{1}{f\_S}\right)\right)} P\left(\mathbf{r} - \frac{1}{f\_S}\right) \right). \tag{4.20}$$

#### **4.2.4.4 Number of Mode Changes and Lifetime**

Since an ESS is not fully charged and discharged with every cycle, the number of times that each ESS changes from charging to discharging mode in each day is calculated in

this work, rather than calculating the number of full cycles. The number of mode changes is simply the number of zero crossings of allocated power profile the ESS (()).

Since the lifetime and the number of cycles is limited for lithium-ion batteries, several contributing factors to the aging this type of ESS are considered, including [17], [199]:


### **4.2.5 Allocating Power Profiles**

The use of a low-pass filter to separate high- and low-frequency components of the variations in the active power has become a standard approach for a simultaneous sizing of several storage technologies. This is particularly the case for sizing a batterysupercapacitor hybrid ESS, where the supercapacitor can increase the power capability of the battery, and preserve its lifetime [213], [214]. The choice of the cut-off frequency of the filter can significantly impact the sizing outcome [197], [214], which can be any value in the range of (0, s⁄2), in which <sup>s</sup> is the sampling frequency of the measured data. Typical discharge times of each storage technology [215], market intervals [216], and optimization frameworks [192] have been used for selecting the filter's cut-off frequency. However, there is no generalized approach for selecting the cut-off frequency, which can be applied to all applications. In our use case, a relatively low cut-off frequency allocates most high-frequency components to the type 2 ESS, but also increases its nominal capacity, which can be quite costly considering the low energy density of these technologies. With a relatively high cut-off frequency, the type 1 ESS has to cover many fast variations in the power profile, which limits the advantages gained by using the second ESS. In this work, a simple optimization framework is introduced for selecting the cut-off frequency, which aims to simultaneously minimize the capacity of the high power density ESS, while reducing the maximum ramp rate and power of the high energy density ESS. This optimization is based on the following assumptions:


With the above-mentioned assumptions, a simple objective function is defined as

$$\text{minimize}\qquadObj = \mathbf{k}\_1 \frac{\frac{E\_{\text{n2}}}{E\_{\text{n1}}}}{\max(\frac{E\_{\text{n2}}}{E\_{\text{n1}}})} + \mathbf{k}\_2 \frac{\frac{P\_{\text{n1}}}{P\_{\text{n2}}}}{\max(\frac{P\_{\text{n1}}}{P\_{\text{n2}}})} + \mathbf{k}\_3 \frac{\frac{R\_{\text{n1}}}{R\_{\text{n2}}}}{\max(\frac{R\_{\text{n1}}}{R\_{\text{n2}}})}.\tag{4.21}$$

In Eq. (4.21), the nominal capacity, the nominal power, and the maximum ramp rate are calculated according to Eq. (4.18)-(4.20) for each type of ESS, which are denoted with 1 and 2 for the type 1 and type 2 ESS, respectively. Also, k<sup>1</sup> , k<sup>2</sup> , and k<sup>3</sup> are weighting coefficients, which are assumed to be equal to one in this work, but can be adjusted when required. The objective function of Eq. (4.21) aims to minimize simultaneously the ratio of capacity of the type 2 the ESS to the capacity of the type 1 ESS, the ratio of the power of the type 1 ESS to the power of the type 2 ESS, and the ratio of the maximum ramp rate of the type 1 ESS to the one of the type 2 ESS. All these three factors are normalized using their maximum value within the range (0, s⁄2). To find the optimum cut-off frequency, a simple gradient search method is used, similar to [197], in which the cut-off frequency is increased from its minimum value up to the point that is no longer decreasing. The frequency, at which the objective function is at its minimum is used as the cut-off frequency of the filter for decoupling the standard patterns. These values are listed in Table 4.2 for the standard pattern of each bus using the 1-second and 1-minute data resolutions. As the discharge time of the type 2 ESS is often much shorter than 10 minutes, their characteristics is not calculated using the 10-minute data.

The optimum cut-off frequencies shown in Table 4.2 are used to decompose the standard patterns into high- and low-frequency components and to allocate these power profiles to each ESS type. Using this approach, the nominal capacity of type 2 ESS is minimized, while simultaneously reducing the power and the maximum ramp rate of type 1 ESS for a prolonged lifetime. Figure 4.18 depicts the change in the objection function of Eq. (4.21) for each bus with the increase in the cut-off frequency of the filter.

As discussed, using the low-pass filter, the high-frequency components of the active power variations are covered by a type 2 ESS. As an example, Figure 4.19 shows the decomposed power profile used for sizing each type of ESS at bus 1. As seen, the shortterm power peaks, for instance, the 40-kW short pulsed power after 8 a.m., are covered by the type 2 ESS, while the type 1 ESS covers only the slower variations in power.



Figure 4.18: Changes in the objective function with the increase of the cut-off frequency of the filter at each bus using the data with 1-second resolution.

Figure 4.19: The allocated power to each type of ESS at bus 1 by applying a low-pass filter with the optimum cut-off frequency.

#### **4.2.6 Sizing Outcome**

Using the standard patterns of each bus and the allocated power profile to each ESS by applying the low-pass filter with the optimum cut-off frequency, the ESS characteristics has been calculated for each bus and each data resolution from Eq. (4.18)-(4.20). The results are shown in detail in Table 4.3.


Table 4.3: Sizing results according to the standard patterns at each bus for the three different data resolutions.

Considering the sizing outcome for bus 1 using the 1-second data as an example, the rated capacity and power for the type 1 ESS is calculated to be 615.12 kWh and 94.87 kW, respectively. For the type 2 ESS with a high power density (e.g., a FESS), the nominal capacity and power are estimated to be 2 kWh and 39 kW. More importantly, using the type 2 ESS, the maximum ramp rate of type 1 ESS is calculated to by only 4.84 kW/s, while the type 2 ESS has to provide a ramp rate of around 35 kW/s. Without the type 2 ESS, this high ramp rate has to be covered by the type 1 ESS. Moreover, the number of mode changes, where the ESS changes from charging to discharging mode is significantly reduced for the type 1 ESS, as most of the small variations in power are covered by the type 2 ESS. Similar results are obtained for other buses, as seen in Table 4.3. The maximum ramp rate and the number of mode changes for the four buses using the 1-second data are also plotted in Figure 4.20 for an easier comparison.

Figure 4.20: (a) The maximum ramp rate, and (b) the number of mode changes per day for both types of ESS using the 1-second data.

For buses 2-4, a type 2 ESS with capacity of less than 1 kWh is required for smoothing out the power profile of type 1 ESS. When comparing the results of the four buses, it can be noted that in buses 3 and 4, which are residential areas with a high share of PV generation, the type 2 ESS requires less ramping rate and a smaller number of mode changes in comparison to bus 1, where industrial loads are connected. Therefore, it can be seen that the variability of industrial loads such as the ones connected at bus 1 can be more demanding for an ESS, in comparison to the variability of PV generated power.

In the following the effect of data resolution, the cut-off frequency of the low-pass filter, and selecting patterns other than the most reoccurring one on the sizing outcome is investigated.

#### **4.2.6.1 The Effect of Data Resolution on the Sizing**

As shown in Table 4.3, the ESS characteristics has been calculated using the data with 1 second, 1-minute, and 10-minute resolutions. It can be observed that the nominal power and capacity of a type 1 ESS are more or less in the same range in most cases, when lower data resolutions are used. The exception is using the 10-minute data for bus 2, since the pattern detected using the 10-minute data is quite different from the pattern detected using other data resolutions. Therefore, it can be concluded that low resolution data up to 10-minutes does not necessarily lead to large errors in estimating the capacity and power of type 1 ESS such as a lithium-ion battery, for the application of peak shaving and balancing the load and PV generation. However, the maximum ramp rate and the number of times the ESS changes from charging to discharging mode and vice versa, cannot be accurately estimated when low data resolutions are used. Hence, using high resolution data, such as 1-second data, should be preferred, when available.

For sizing the type 2 ESS for power smoothing applications, it is necessary to use data with the highest resolution. Even using 1-minute data can lead to extreme over-sizing of the ESS, as in the case of buses 3 and 4, where the estimated capacity using the 1-minute data is twice the estimated value using the 1-second data. Similar to type 1 ESS, using low resolution data leads to a large underestimation of the required ramp rate and number of mode changes per day for the type 2 ESS.

#### **4.2.6.2 The Effect of Filter's Cut-off Frequency**

The sizing result of Table 4.3 were based on the optimum cut-off frequency calculated from the objective function of Eq. (4.21), and given in Table 4.2. Figure 4.21 shows how the characteristics of the type 1 and 2 ESS alters for bus 2, as an example, if other values are used for the cut-off frequency of the low-pass filter.

Figure 4.21: The variations in the characteristics of the type 1 and type 2 ESS with respect to the change in the cut-off frequency of the low-pass filter for bus 2. The optimum frequency is the red dot.

As seen from Figure 4.21, with an increase in the filter's cut-off frequency, the power and the maximum ramp rate of the type 1 ESS increases, while its capacity remains relatively unchanged. On the contrary, the capacity and power of type 2 ESS decreases drastically with the raising cut-off frequency, as less power is allocated to this storage technology, up to reaching the Nyquist frequency, where all frequency components are allocated to the type 1 ESS. The optimum cut-off frequency, shown with the red marker in Figure 4.21, is the point where the capacity of the type 2 ESS is relatively low (below 1 kWh), while simultaneously having a relatively low ramp rate and power demanded from the type 1 ESS.

#### **4.2.6.3 The Effect of Selecting Different Patterns**

As shown in Figure 4.17, the motif discovery algorithm can possibly generate several different patterns, which are also reoccurring, but with lower number of reoccurrences. Table 4.4 shows the sizing outcome for bus 2 using the 1-minute data, if other patterns are selected. As discussed before and shown in Figure 4.17, pattern 2 is the pattern of a weekend, while pattern 3 corresponds to the first day after a weekend or holiday.


Table 4.4: Sizing results using three different standard at bus 2 using the 1-minute data.

It can be perceived from Table 4.4 that selecting different patterns can lead to different sizing outcomes for the type 1 ESS. When using pattern 2, which is reoccurring pattern during the weekends or holidays, the nominal power of the type 1 ESS is largely underestimated. When pattern 3 is used, which is the pattern of the first working day of the week, the calculated nominal capacity is significantly higher, due to the different power requirement during the first 6 hours of the day (see Figure 4.17). This shows how a random and arbitrary choice of the day for the sizing study, which indeed could be a weekend or a first day of the week, can lead to very different results, in particular for sizing the type 1 ESS. Therefore, such an approach should be avoided.

For the type 2 ESS, the sizing outcomes are not considerably different, when other patterns are used. This is due the fact that the high-frequency components above the filter's cut-off frequency have been more or less consistent in all patterns in our use case.

### **4.2.7 Evaluation of the Proposed Sizing Method**

In order to evaluate the proposed sizing approach using standard patterns, the ESS with the characteristics derived from the standard patterns are tested throughout the whole data set. Since not all days follow exactly the same pattern, and the 80%-quantile of the most reoccurring pattern is selected for the storage sizing, it is evident that the ESS with the characteristics based on the standard pattern will not be able to fulfil its purposes at all times and during all days. However, since the ESS is sized based on the daily pattern that has been repeated the most, it should be capable to fulfil its purposes for most of the days within the data set. The performance of the ESS is evaluated based on how well the peak shaving and power smoothing goals can be fulfilled using the two types of ESS throughout the whole measurement period.

In order to see the effect of using type 2 ESS, such as a FESS, for power smoothing applications, the results are shown for the 1-second data only. Also, bus 3 is selected for demonstrating the performance of the storage systems, as it represents a typical duck curve with a high share of PV generation. Similar results can be shown for the other buses and using other data resolutions. Figure 4.22 shows the active power drawn from the medium voltage grid at bus 3 during the period of two weeks, which is the length of the original data set used for detecting patterns.

The results are shown for three different cases:


Within the two weeks of measurement data, the capacity, nominal power, and the maximum ramp rate calculated from the standard patterns are used as the constraints for charging and discharging of both types of ESS. Also, the state of charge of both types of storage systems are maintained within the range of the state of charge discussed in subsection 4.2.4, which is also illustrated in Figure 4.23. A simple on/off controller is used for the peak shaving application using the type 1 ESS, which stores the PV generation starting from 10:00 until the ESS is fully charged, and discharges it from 18:00, until the storage is empty. These times are also selected based on the detected standard pattern at this bus. For the type 2 ESS, a continuous operation is considered in order to continuously filter out short-term power variations.

Figure 4.22: The active power drawn from the medium voltage grid at bus 3, with and without the type 1 ESS, and in combination with a type 2 ESS, such as a FESS, (a) during the whole two weeks, (b) during the 5th day, and (c) between 22:00 and 23:00 of the 5th day.

Figure 4.23: The state of charge variations during the two weeks for both ESS types, with the characteristics derived from the standard patterns.

From Figure 4.22, it is evident that:


It can be concluded, that each ESS with the characteristics derived from the standard patterns can fulfil its purposes for most of the days within the data set, showing the effectiveness of the proposed method.

# **4.3 Summary**

In this chapter, two novel data-driven solutions were introduced for the problem of allocation and sizing of an ESS.

The first section presented an allocation algorithm for a FESS for power smoothing applications using the concept of mutual information. A data-driven method based on the concept of mutual information was proposed for estimating the relative voltage sensitivity to active and reactive power changes. Unlike the classical linearized approach, the proposed approach does not require a grid model, and can capture system nonlinearities and the changes in the system operation point. The suggested method was applied for the allocation of a FESS in a real grid in southern Germany. By collecting and using measured data at the most attractive candidates, it has been shown that the proposed approach can successfully find the locations with a high voltage sensitivity to active changes. These locations can be suitable for the installation of a FESS, if sharp active power variations, for example, caused by pulsed power loads, are also present.

In the second section, a new data-driven method is introduced for sizing different types of ESS, including a high power density ESS, such as a FESS, for power smoothing applications and a high energy density ESS, such as a lithium-ion battery, for reducing the evening peak and avoiding reverse power flow. The main contribution of the proposed method is finding reoccurring daily consumption patterns using the motif discovery algorithm, and using it for the sizing calculations, rather than using an arbitrary day, or reducing the data resolution for the sizing study, as commonly found in literature. Motif discovery can find similar, and yet, not the same patterns using the SAX representation of the data, and generate a representative day for the time series, which can be used for the storage sizing. The detected patterns using motif discovery were used for deriving the characteristic of two types of ESS. Using a simple optimization framework, an optimum cut-off frequency is selected for applying a low-pass filter on the patterns in order to decompose them into low- and high-frequency components, and allocating them to different types of ESS. The nominal capacity, nominal power, and the maximum ramp rate is calculated for each type of ESS using the standard patterns, detected from the measurement data collected at four different low voltage distribution grids in southern Germany with three different data resolutions. It has been shown that the ESS with characteristic derived from only the standard patterns can effectively reduce the evening peak, the reverse power flow, and short-term active power fluctuations for most of the days throughout the whole data set. Moreover, it has been shown, that high-resolution data should be preferred to accurately calculate all the characteristics of an ESS, in particular for power smoothing applications.

In summary, this chapter presented two examples on how techniques and solutions from computer science and information theory can be applied for solving practical problems in power systems. With the increasing number of measurement devices in particular at the low voltage distribution grids, data-driven techniques can become a major tool for providing novel solutions to the challenges in design and operation of low voltage grids.

# **5 Modelling of High-Speed Flywheel Energy Storage Systems (FESS) for Real-time Simulations**

In this chapter, the model of a high-speed Flywheel Energy Storage System (FESS) is presented and simulated in real-time, which considers the losses and auxiliary power requirements of the FESS. The focus of the modelling is on the behavior of the FESS for frequency and voltage support at the low voltage level. The model of the FESS is simulated in real-time with a simulation step of 10 µs, and the results are compared with offline simulation results. The developed real-time simulation model can be used for controller hardware-in-the-loop testing of new controllers, or for the preliminary studies required for Power Hardware-in-the-Loop (PHIL) testing of the FESS, as discussed in Chapter 3. The model of the FESS in this work is based on the structure and design of a specific commercial 60 kW high-speed FESS. However, the model can be easily adjusted for any other systems.

# **5.1 Potential Applications of Real-time Simulation of a FESS**

As discussed in Chapter 3, real-time simulations of new power systems components can accelerate their design, development, and performance verification [96], [100]. A component model, running in real-time, can serve as a *Digital Twin* [99] of that component, which can be used to study its interaction with other hardware. Therefore, engineers have developed models of power systems components including wind turbines [111], [217], photovoltaic systems [218], and superconducting fault current limiters [219], [220] for real-time simulations.

Moreover, a real-time simulation of a power system component together with a power amplifier can be used to build a full-scale emulator of the component. These emulators are often the safer and the cheaper alternatives to the real hardware, for instance, in the case of hazardous lithium-ion batteries or the costly high-speed FESS. Lithium-ion batteries have been emulated at the cell level [221], as small battery energy storage system for household applications [222], and as grid-scale systems for grid-wide applications [223]. Today, commercial and non-commercial emulators are available for many power system components such as wind turbines [143], [224], PV systems [225], [226], and electrical machines [227], which are used for various applications, including control and converter testing [173], [228], [229].

In addition, as discussed in Chapter 3, real-time simulation of the hardware under test is a recommended step, prior to a PHIL testing. Although offline simulations can also provide useful insights [118], there are advantages of having the model of the hardware under test running in real-time. This includes simulating the hardware model in the same grid model, in which the behavior of the real hardware is later evaluated using PHIL testing.

To our knowledge, the real-time implementation of a FESS model has not yet been carried out before. In this chapter, a model of a high-speed FESS is developed and simulated in real-time, and the results are compared against offline simulation results. In the following chapter, this model is validated using the experimental results obtained from the PHIL testing of a commercial FESS.

# **5.2 Modeling of a High-speed Flywheel Energy Storage System (FESS)<sup>1</sup>**

The model of a FESS comprises the model of each of its individual components. Figure 5.1 illustrates the structure of a high-speed FESS at a component level. The structure of the high-speed FESS presented in Figure 5.1 is based on the design of a 60 kW/3.6 kWh commercial high-speed FESS. The inside view of the container holding the 60 kW FESS and its required auxiliary components is shown in Figure 5.2.

As seen in Figures 5.1 and 5.2 and discussed in Chapter 2, a high-speed FESS consists of a Permanent Magnet Synchronous Machine (PMSM) coupled to a high-inertia rotor, two Voltage Source Converters (VSC) as the Machine-side Converter (MSC) and Grid-side Converter (GSC) with their corresponding controllers, and filters on the AC side of each converter. Due to the extremely low inductance values of the PMSM used in this highspeed FESS, it is necessary to have an LC filter between the PMSM and the machine-side converter to reduce current and torque ripples [230].

<sup>1</sup> Parts of this chapter have been presented in the following publications:

<sup>1)</sup> S. Karrari, M. Noe, and J. Geisbuesch, "High-speed Flywheel Energy Storage System (FESS) for Voltage and Frequency Support in Low Voltage Distribution Networks," in 3rd IEEE International Conference on Intelligent Energy and Power Systems (IEPS 2018), Kharkov, 2018, pp. 176–182.

<sup>2)</sup> S. Karrari, M. Noe, and J. Geisbuesch, "Real-time simulation of high-speed Flywheel Energy Storage System (FESS) for distribution networks," in the 9th ACM International Conference on Future Energy Systems (ACM e-Energy 2018), Karlsruhe, 2018.

Figure 5.1: Structure of a high-speed Flywheel Energy Storage System (FESS).

Figure 5.2: The inside view of the container of the 60 kW high-speed FESS.

In principle, the model of a FESS is similar to the model of a variable-speed wind turbine equipped with a PMSM [107]. However, the FESS has faster dynamics, no mechanical input torque, and the ability to absorb energy. In the following section, the model of each individual component of the FESS is discussed in detail.

## **5.2.1 Modeling of a Permanent Magnet Synchronous Machine**

As mentioned in Chapter 2, Permanent Magnet Synchronous Machines (PMSM) are the most common choice for a high-speed FESS, while a low-speed FESS often uses an asynchronous machine. Due to the absence of the field windings, PMSM are more suitable to be used in the vacuum enclosure of high-speed FESS, where convection cooling of the rotor is not possible. They are also simple and robust in their structure, which leads to a higher reliability. Moreover, these machines have high power-to-weight and torque-to-mass ratios and an easier control design, when compared to several other types of electrical machines [73].

Depending on the arrangement of the permanent magnets in the rotor, there are several types of PMSM, with rather diverse characteristics. Four major types of PMSM are shown in Figure 5.3. In a high-speed FESS, a surface mounted PMSM, which is shown in Figure 5.3(a), is often used, in which the permanent magnets are mounted on the rotor surface in a symmetrical arrangement. This leads to a relatively small stator inductance and a negligible reluctance torque. The other common type is the interior PMSM, shown in Figure 5.3(c) and (d), in which the rotor salience results in a non-uniform air gap flux distribution and a high reluctance torque. A comparison between a surface mounted and an interior PMSM is given in Table 5.1. The interior PMSM and the surface mounted PMSM have different characteristics, however, the mathematical description is the same in both cases [231].


Table 5.1: A comparison between the surface mounted and interior PMSM [231].

For modeling and description of the notations of a PMSM, a simplified schematic of a two-pole surface mounted PMSM is shown in Figure 5.4. Clearly, the modeling and the equations are also valid for a PMSM with more than two poles.

Figure 5.3: Arrangement of the permanent magnets in (a) surface mounted PMSM, (b) surface inset PMSM, (c) interior (I-shaped) PMSM, and (d) interior circumferential PMSM.

Figure 5.4: (a) The equivalent circuit of a PMSM in the three-phase stationary reference frame. (b) Schematic of a two-pole surface mounted PMSM in the dq-rotating reference frame.

The model of a PMSM has been frequently described in the literature, including in [73], [231]–[233]. In this work, a very detailed model is not necessary, as the goal is to represent the behavior of the FESS at the system level for grid applications. Therefore, the following simplifying assumptions have been made for modeling the PMSM [73]:


Moreover, the iron saturation is also neglected in this study due to unavailability of the relevant data from the manufacturer. However, for machines with a high power density that are often used, it is recommended to consider the iron saturation [234]. The losses are neglected in the modeling of the PMSM. However, all losses in the FESS including the copper and iron losses in the PMSM are aggregated and considered, as later described in subsection 5.2.5. The permanent magnets are assumed to be a flux linkage source concentrating all their flux along only one axis [231]. Since a surface mounted PMSM is considered, a uniform air gap with no rotor saliency is assumed.

With the aforementioned assumptions, the PMSM can be represented by the equivalent circuit shown in Figure 5.4(a). The relationship between the voltages and currents of each phase for a balanced system in the three-phase stationary reference frame is governed by [73]

$$
\begin{bmatrix} u\_{\mathbf{a}} \\ u\_{\mathbf{b}} \\ u\_{\mathbf{c}} \end{bmatrix} = \mathbf{R}\_{\mathbf{s}} \begin{bmatrix} l\_{\mathbf{a}} \\ l\_{\mathbf{b}} \\ l\_{\mathbf{c}} \end{bmatrix} + \frac{d}{dt} \begin{bmatrix} \psi\_{\mathbf{a}} \\ \psi\_{\mathbf{b}} \\ \psi\_{\mathbf{c}} \end{bmatrix}, \tag{5.1}
$$

in which,

$$
\begin{bmatrix}
\psi\_{\rm a} \\
\psi\_{\rm b} \\
\psi\_{\rm c}
\end{bmatrix} = \begin{bmatrix}
\mathcal{L}\_{\rm aa} & \mathcal{L}\_{\rm ab} & \mathcal{L}\_{\rm ac} \\
\mathcal{L}\_{\rm ba} & \mathcal{L}\_{\rm bb} & \mathcal{L}\_{\rm bc} \\
\mathcal{L}\_{\rm ca} & \mathcal{L}\_{\rm cb} & \mathcal{L}\_{\rm cc}
\end{bmatrix} \begin{bmatrix}
l\_{\rm a} \\
l\_{\rm b} \\
l\_{\rm c}
\end{bmatrix} + \psi\_{\rm PM} \begin{bmatrix}
\sin\theta\_{\rm m} \\
\sin(\theta\_{\rm m} - \frac{2\pi}{3}) \\
\sin(\theta\_{\rm m} + \frac{2\pi}{3}) \\
\end{bmatrix}.
\tag{5.2}
$$

In Eq. (5.1) and (5.2), <sup>k</sup> is the stator phase voltage (k ∈ {a, b, c}), <sup>k</sup> is the stator phase current, <sup>k</sup> is the stator flux linkage, PM is the permanent magnet flux linkage, Lkk are the self- and mutual inductances of each stator phase, R<sup>s</sup> is the stator phase resistance, and <sup>m</sup> is the mechanical angle of the rotor. According to Eq. (5.1), the voltage of the three-phase windings in a PMSM equals to the sum of the resistive voltage drops and the time derivative of the flux linkages. However, the fluxes and inductances in Eq. (5.2) change nonlinearly with rotation of the rotor. To avoid this, it is common practice to replace the three-phase stationary windings of the stator with two imaginary perpendicular windings, which are rotating with the electrical speed of the rotor. This frame is known as the dq-rotating reference frame and is depicted in Figure 5.4(b). The well-known Park transformation matrix (dq), as shown in Eq. (5.3), is used for mapping the values from the three-phase stationary frame to the dq-rotating reference frame. In Eq. (5.3), <sup>d</sup> and <sup>q</sup> are the d- and q-axis equivalent stator currents, and <sup>0</sup> , which is the zero component current, is assumed to be zero, as a balanced system is considered. As shown in Figure 5.4(b), it is assumed that the d-axis is aligned with the rotating flux produced by the permanent magnets and the current of the first phase (phase a). The qaxis leads the d-axis by 90 degrees in the counterclockwise direction.

$$
\begin{bmatrix}
\dot{l}\_{\text{d}} \\
\dot{l}\_{\text{q}} \\
\dot{l}\_{\text{0}}
\end{bmatrix} = \frac{2}{3} \begin{bmatrix}
\cos\left(\omega\_{\text{e}}t\right) & \cos\left(\omega\_{\text{e}}t - \frac{2\pi}{3}\right) & \cos\left(\omega\_{\text{e}}t + \frac{2\pi}{3}\right) \\
\frac{1}{2} & \frac{1}{2} & \frac{1}{2} \\
\frac{1}{2} & \frac{1}{2} & \frac{1}{2}
\end{bmatrix} \begin{bmatrix}
\dot{l}\_{\text{a}} \\
\dot{l}\_{\text{b}} \\
\dot{l}\_{\text{c}}
\end{bmatrix}.
$$

Since the reference frame rotates with the same frequency as the supplied voltage in steady-state conditions, this frame is also referred to as the synchronous reference frame in the literature [235]. Using the dq-rotating reference frame, sinusoidal variables become DC, which are much easier to analyze and control. In the dq-rotating reference, the PMSM is modeled as [236]

$$
\begin{bmatrix} u\_{\rm d} \\ u\_{\rm q} \end{bmatrix} = \mathcal{R}\_{\rm s} \begin{bmatrix} l\_{\rm d} \\ l\_{\rm q} \end{bmatrix} + \frac{d}{dt} \begin{bmatrix} \psi\_{\rm d} \\ \psi\_{\rm q} \end{bmatrix} + \omega\_{\rm e} \begin{bmatrix} -\psi\_{\rm q} \\ \psi\_{\rm d} \end{bmatrix}, \tag{5.4}
$$

in which,

$$
\begin{bmatrix}
\psi\_{\rm d} \\
\psi\_{\rm q}
\end{bmatrix} = \begin{bmatrix}
\mathcal{L}\_{\rm d} \dot{\iota}\_{\rm d} + \psi\_{\rm PM} \\
\mathcal{L}\_{\rm q} \dot{\iota}\_{\rm q}
\end{bmatrix} \tag{5.5}
$$

and <sup>d</sup> and <sup>q</sup> are the d- and q-axis stator voltage, L<sup>d</sup> and L<sup>q</sup> are d- and q-axis stator inductance, <sup>d</sup> and <sup>q</sup> are the stator flux in the d- and q-axis directions, and <sup>e</sup> is the rotor electrical speed. As seen from Eq. (5.5), using this approach, the flux values are no longer a function of the rotor position. As saturation is neglected, a linear relationship is assumed between the d- and q-axis flux linkages and the stator currents, as shown in Eq. (5.5). A detailed modeling of the PMSM considering saturation, nonlinear cross coupling, and angular dependencies can be found in [227].

Equations (5.4) and (5.5) can be schematically illustrated by two equivalent circuits, one for each of the d- and q-axis directions, as shown in Figure 5.5. In this figure, the voltage sources represent the induced voltages caused by the rotation of the machine.

Figure 5.5: Equivalent circuits of a PMSM in the (a) d-axis and the (b) q-axis.

For calculating the speed variations of the rotor, the mechanical part of the PMSM is modeled using the Eq. (5.6)-(5.10), considering the fact that in a FESS no mechanical torque input is applied to the rotor.

$$
\omega\_{\mathbf{e}} = n\_{\mathbf{p}} \omega\_{\mathbf{m}}.\tag{5.6}
$$

$$
\omega\_{\rm m} = \frac{d\theta\_{\rm m}}{dt},
\tag{5.7}
$$

$$\begin{aligned} \mathbf{g} &= \frac{3}{2} n\_{\mathbf{p}} \left( \boldsymbol{\psi}\_{\mathbf{d}} \boldsymbol{i}\_{\mathbf{q}} - \boldsymbol{\psi}\_{\mathbf{q}} \boldsymbol{i}\_{\mathbf{d}} \right) = \overbrace{\frac{3}{2} n\_{\mathbf{p}} \boldsymbol{i}\_{\mathbf{d}} \boldsymbol{i}\_{\mathbf{q}} \left( \mathbf{L}\_{\mathbf{d}} - \mathbf{L}\_{\mathbf{q}} \right)}^{\text{Reluctance torque}} + \overbrace{\frac{3}{2} n\_{\mathbf{p}} \boldsymbol{\psi}\_{\mathbf{p}} \boldsymbol{\psi}\_{\mathbf{p}} \boldsymbol{i}\_{\mathbf{p}}}^{\text{synfations torque}}, \end{aligned} \tag{5.8}$$

$$\begin{aligned} \tau\_{\mathbf{e}} &= \frac{1}{2} n\_{\mathbf{p}} (\psi\_{\mathbf{d}} i\_{\mathbf{q}} - \psi\_{\mathbf{q}} i\_{\mathbf{d}}) = \frac{1}{2} n\_{\mathbf{p}} i\_{\mathbf{d}} i\_{\mathbf{q}} (\mathcal{L}\_{\mathbf{d}} - \mathcal{L}\_{\mathbf{q}}) + \frac{1}{2} n\_{\mathbf{p}} \psi\_{\mathbf{p} \mathbf{M}} i\_{\mathbf{q}} \\ \mathcal{J}\_{\mathbf{f}} \frac{d\omega\_{\mathbf{m}}}{dt} &= \tau\_{\mathbf{e}} - \mathcal{D}\_{\mathbf{f}} \omega\_{\mathbf{m}} \end{aligned} \tag{5.9}$$

$$\begin{aligned} \mathbf{J}\_{\mathbf{f}} \frac{\mathbf{\cdot} \cdot \mathbf{m}}{\mathbf{d}t} &= \mathbf{\tau}\_{\mathbf{e}} - \mathbf{D}\_{\mathbf{f}} \boldsymbol{\omega}\_{\mathbf{m}}, \\ E\_{\mathbf{F}\mathbf{m}} &= \frac{1}{\mathbf{I} \cdot \boldsymbol{\omega}\_{\mathbf{m}}^2} \mathbf{I} \end{aligned} \tag{5.10}$$

$$E\_{\rm FW} = \frac{1}{2} \mathbf{J}\_{\rm f} \boldsymbol{\omega}\_{\rm m}^{2}.\tag{5.10}$$

In these equations, <sup>m</sup> is the mechanical angular velocity of the rotor, <sup>e</sup> is the electrical torque, <sup>p</sup> is the number of pole pairs, D<sup>f</sup> is the friction coefficient, FW is the energy stored in the rotor, and J<sup>f</sup> is the combined moment inertia of the PMSM and the flywheel. Equation (5.9) determines the variations in the speed of the rotor or flywheel, from which the changes in its energy content can be calculated according to Eq. (5.10). The term Df<sup>m</sup> in Eq. (5.9) is referred to as the drag toque, which is minimized in a highspeed FESS with the use of a vacuum enclosure and magnetic bearings. Equation (5.10) resembles the swing equation of conventional synchronous generators. This similarity is used later in Chapter 7 to introduce an inertia emulation controller for a FESS in lowinertia grids.

The electrical torque in Eq. (5.8), consists of two parts; the reluctance torque and the synchronous torque. In a surface mounted PMSM, the d- and q-axis stator inductances (L<sup>d</sup> and Lq) are almost equal, making the reluctance torque nearly zero. Therefore, Eq. (5.11) can be further simplified to

$$
\pi\_{\mathbf{e}} = \frac{3}{2} \mathbf{n}\_{\mathbf{p}} \psi\_{\mathbf{PM}} \mathbf{i}\_{\mathbf{q}}.\tag{5.11}
$$

Since the permanent magnet flux is assumed to be constant<sup>2</sup> , it can be seen from Eq. (5.11) that the electromagnetic torque can be controlled through the q-axis component of the stator current. This simplicity cannot be achieved when using other reference frames such as the stationary αβ-frame.

The active (<sup>e</sup> ) and reactive power () of a PMSM in the dq-rotating reference frame is calculated using Eq. (5.12) and (5.13), which is also known as the instantaneous power theory [237].

$$P\_\mathbf{e} = \frac{3}{2} \{ u\_\mathbf{d} i\_\mathbf{d} + u\_\mathbf{q} i\_\mathbf{q} \}. \tag{5.12}$$

$$Q = \frac{3}{2} \{-u\_d i\_q + u\_q i\_d\}.\tag{5.13}$$

The model of the PMSM in the dq-rotating reference frame can also be rewritten as a nonlinear state-space model by defining ̇d, ̇q, and <sup>e</sup> as the state variables, and and as the system input, which result

$$\frac{d}{dt}\begin{bmatrix}\frac{l\_{\text{L}}}{\mathbf{L}\_{\text{d}}} & 0 & 0\\ \frac{-\mathbf{R}\_{\text{s}}}{\mathbf{L}\_{\text{d}}} & \frac{-\mathbf{R}\_{\text{s}}}{\mathbf{L}\_{\text{q}}} & \frac{-\boldsymbol{\psi}\_{\text{PM}}}{\mathbf{L}\_{\text{q}}}\\ 0 & \frac{3}{\mathbf{L}\_{\text{q}}} & \frac{-\boldsymbol{\psi}\_{\text{PM}}}{\mathbf{L}\_{\text{q}}}\\ 0 & \frac{3}{2}\frac{\mathbf{n}\_{\text{p}}^{2}}{\mathbf{l}\_{\text{f}}}\boldsymbol{\psi}\_{\text{PM}} & \frac{-D}{\mathbf{l}\_{\text{f}}}\end{bmatrix}\begin{bmatrix}\frac{1}{\mathbf{L}\_{\text{d}}}\\ \frac{1}{\mathbf{L}\_{\text{d}}}\\ \frac{-\mathbf{L}\_{\text{q}}}{\mathbf{L}\_{\text{d}}}\end{bmatrix} + \begin{cases} \frac{\mathbf{L}\_{\text{q}}}{\mathbf{L}\_{\text{d}}}\boldsymbol{\omega}\_{\text{d}}\boldsymbol{i}\_{\text{q}}\\ \frac{-\mathbf{L}\_{\text{q}}}{\mathbf{L}\_{\text{d}}}\boldsymbol{\omega}\_{\text{d}}\boldsymbol{i}\_{\text{d}}\\ \frac{-\mathbf{R}\_{\text{q}}}{\mathbf{L}\_{\text{d}}}\boldsymbol{\omega}\_{\text{d}}\boldsymbol{i}\_{\text{d}}\\ \frac{3}{2}\frac{\mathbf{n}\_{\text{p}}}{\mathbf{L}\_{\text{f}}}\boldsymbol{i}\_{\text{d}}\boldsymbol{i}\_{\text{d}}\left(\mathbf{L}\_{\text{d}}-\mathbf{L}\_{\text{q}}\right)\end{bmatrix}.\tag{15.14}$$

Knowing the exact values of the machine parameters is mandatory for an accurate simulation of its dynamic and steady-state behavior. The permanent magnet flux (PM)

<sup>2</sup> Minor variations due to temperature changes may be present in practice, which are neglected in this work.

can be easily calculated using an open-circuit test with a reasonable accuracy. However, determining the precise values for d- and q-axis stator inductances (L<sup>d</sup> and Lq) can be more complex, as saturation and cross coupling effects are also affecting the results. Various analytical, finite element, and experimental methods can be used to estimate these values, as described fully in [238], [239]. In this work, these values are provided by the manufacturer and are given in Table 5.2.


Table 5.2: Parameters of the PMSM of the high-speed FESS according to the manufacturer.

### **5.2.2 Modeling of Voltage Source Converters**

As shown in Figure 5.1, a FESS is controlled via two voltage source converters in a backto-back configuration. As discussed in Chapter 2, several converter topologies can be used in a FESS, but the two-level and three-level converters are the most common topologies found in commercial systems. In the high-speed FESS in this work, both converters have a two-level topology.

Depending on the intended study, there are various well-known converters models that can be selected, which differ mainly on how detailed the switches and semiconductors are represented. Comparative studies of different converter models for various applications can be found in [110], [240], [241]. Clearly, more detailed models have a higher accuracy, but at the cost of a higher computational effort. A small computation time is of crucial importance for real-time simulations, in which all calculations are required to be completed within a short simulation time step. Therefore, choosing the appropriate level of details is necessary to limit the model computation time, while simultaneously having an acceptable model accuracy. Average model, switching model, and detailed models are among the most commonly used converter models. These well-known models are explained briefly in this section, followed by a discussion on choosing the right model for the real-time simulations in this work. It should be noted that there are several other lessknown converter modeling techniques, such as the sub-cycle average models [242], which are not discussed here.

#### **5.2.2.1 Average Model**

The average model is based on the assumption that the switching frequency of the converter is much greater than the fundamental frequency of the grid. In this model, the switching transitions, harmonics, and ripples are removed from the voltages and currents by averaging over one switching period. Nevertheless, lower frequency components of voltages and currents are preserved.

In the average model, all switches are assumed to be operated in the continuous conduction mode [235]. The semiconductors like IGBTs and diodes are not explicitly modeled, and the converter on the AC side is simply modeled as a three-phase controlled AC voltage source with a small internal resistance using [243]

$$
\mu\_{\rm c\downarrow} = m\_{\rm l} \frac{\upsilon\_{\rm DC}}{2} - r\_{\rm on} \iota\_{\rm c\downarrow}.\tag{5.15}
$$

In Eq. (5.15), <sup>j</sup> is the modulation index of each phase (j ∈ {a, b, c}), which is determined by the converter controller, DC is the voltage of the DC side, cj and cj are the phase voltage and current of the AC side, respectively, and on represents the on-state resistance of the switches [243]. It is important to ensure that |<sup>j</sup> | ≤ 1 at all times in order to guarantee a linear Pulsed Width Modulation (PWM). Equation (5.15) can be rewritten in the dq-rotating reference frame, rotating at the fundamental frequency of the AC side voltage using [243]

$$
\begin{bmatrix} \boldsymbol{u}\_{\rm cd} \\ \boldsymbol{u}\_{\rm cq} \end{bmatrix} = \frac{\boldsymbol{U}\_{\rm DC}}{2} \begin{bmatrix} \boldsymbol{m}\_{\rm d} \\ \boldsymbol{m}\_{\rm q} \end{bmatrix} - \boldsymbol{r}\_{\rm on} \begin{bmatrix} \boldsymbol{i}\_{\rm cd} \\ \boldsymbol{i}\_{\rm cq} \end{bmatrix}. \tag{5.16}
$$

In Eq. (5.16), cd and cq are the d- and q-axis voltage of the AC side of the converter, cd and cq are the d- and q-axis current, and <sup>d</sup> and <sup>q</sup> are the d- and q-axis component of the modulation index, respectively.

For modeling the DC side, a lossless converter is assumed (losses are later aggregated in subsection 5.2.5). Therefore, the current of the DC side of the converter is calculated using

$$I\_{\rm DC} = \frac{P\_{\rm AC}}{U\_{\rm DC}} = \frac{\frac{3}{Z} \left( u\_{\rm cd} l\_{\rm cd} + u\_{\rm cq} l\_{\rm cq} \right)}{U\_{\rm DC}}.\tag{5.17}$$

101

According to Eq. (5.17), the DC side of the converter is modeled as a controlled DC current source, in which the current is a function of the active power flowing through the AC side (AC). An equivalent circuit of the average model of a voltage source converter is shown in Figure 5.6.

Figure 5.6: The equivalent circuit of the average model of a voltage source converter.

Since all switching transitions, harmonics, and ripples are neglected, the average model requires a very low computational effort. It is the most common model used for designing converter controllers, as it provides a simple physical insight. Controllers can be included in full detail, considering the inner, e.g., current, and outer, e.g., power or DC-link voltage, control loops. From the grid perspective, since the high-frequency harmonics are mostly attenuated by filters, the average model can be suitable for integration studies of distributed energy resources into power systems. It can be used for transient stability studies of fast phenomena in the order of tens to hundreds of milliseconds [110]. However, this model is not suitable for the investigation of harmonics or internal converter analyses.

#### **5.2.2.2 Switching Model**

In the simplest form of switching model, the switches are modeled as ideal on/off switches. The switching transitions such as dead times are ignored based on the assumption that they take only a small fraction of the switching period. Each switch is controlled by a switching function, which has the value of one, when the switch is on, and zero, when the switch is off. On-state losses can be added as a resistor in series with the ideal switches. Therefore, the voltage of each converter phase using the switched model is calculated using [243]

$$\mathcal{U}\_{\rm cj} = \left(\mathcal{S}\_{\rm mj} - \mathcal{S}\_{\rm nj}\right)\frac{U\_{\rm DC}}{2} - r\_{\rm on}i\_{\rm cj} \tag{5.18}$$

where mj and nj are the switching function of two complementary switches for each converter leg, and

$$\mathcal{S}\_{\text{m}\downarrow} + \mathcal{S}\_{\text{n}\downarrow} = 1\tag{5.19}$$

at all times, in order to avoid short-circuiting of the DC-link. The switching functions are the gate signals from the pulsed width modulation controller. Modeling of the DC side is similar to the one from the average model by calculating the active power from the three phases voltages and current. The equivalent circuit for a voltage source converter using the switching model is shown in Figure 5.7

Figure 5.7: The equivalent circuit of the switching model of a voltage source converter.

In a more detailed approach, switches can be modeled as two-value resistors, having a high resistance when they are in the off-state (in the range MΩ) and a low resistance, when they are in the on-state (in the order of mΩ). The snubber circuit and the antiparallel diodes are also modeled in the same manner. The choice between high and low resistance values for the switches is made on the basis of the gate signals received from the modulation controller, but this choice for the anti-parallel diodes are made based on the instantaneous voltage and current values.

Switching models can represent voltage and current harmonics, which is a significant advantage compared to the average model. However, small simulation time steps are required usually for simulating such models to avoid numerical instability, in particular for high switching frequencies.

#### **5.2.2.3 Detailed Model**

There are different types of detailed converter models in literature, as described in [240]. In one approach, switches are represented by variable resistors, which can either follow a simplified linear function or a nonlinear one according to the voltage-current characteristics of the semiconductors, given in the manufacturer's datasheets [107]. The detailed modeling of converters can be even more complex, in which all the parasitic capacitors of the switches, stray inductances, and the packaging effects are considered.

Simulating such models usually requires extremely small simulation time steps to accurately represent fast switching events and solve the stiff differential equations of the system with high order solvers. The combination of small simulation time steps and high model complexity leads to a long run time for these models. Therefore, such models are often used for a specific part of the system and for a short time period of a few switching cycles. They are mostly used for specific applications, such as simulating internal converter faults, studying switching transitions, switching losses and electromagnetic interference (EMI). This level of details is not required for power system integration studies or designing and testing control systems [110], [240], [244], which is the intention in this work.

#### **5.2.2.4 Selecting the Right Converter Model for Real-time Simulations**

As discussed in Chapter 3, in real-time simulations, equations should be solved within a fixed simulation time step. The average model can easily be simulated in real-time, independent of the switching frequency, as variations within each switching period are neglected. Real-time simulation of a switching model for converters can be challenging for high switching frequencies. As discussed in subsection 3.1.2, in practice, the simulation sampling frequency is usually chosen in the range of 50 to 100 times of the switching frequency for an accurate simulation of converter switching models [111], [245]. Furthermore, in the switching model using two-value resistors, there is a sudden change between a high-value resistance and a low-value one. Hence, an even smaller time step may be required to avoid numerical instabilities, when using the fixed-step noniterative solvers of the real-time simulators [107]. Therefore, real-time simulation of converters with a high PWM carrier frequency using the switching models may require simulation time steps of less than few microseconds, which is not possible using CPUbased real-time simulators. In this case, the model has to be simulated on an FPGA with a submicrosecond simulation time step.

In this work, two voltage source converters with a PWM carrier frequency of 8 kHz have to be simulated in real-time. This requires a simulation step of 1.25-2.5 µs for accurate real-time simulation of the switching models according to [111], [245], which is not possible using CPU-based real-time simulators. However, as suggested in [110], for applications such as frequency and voltage control in power systems, which are the applications intended in this work, the use of average converter models can be adequate. Therefore, average converter models using Eq. (5.15)-(5.18) are used for the modeling of the FESS. This also allows using larger simulation time steps, which means more components can be simulated in real-time in a single model. It should be noted that simulating the switching model of these converters on the FPGA or with the help of switching event compensation algorithms [246] is also possible, but deemed unnecessary for this application.

### **5.2.3 Grid-side Converter Controller**

Figure 5.8 depicts the block diagrams of the grid-side converter and machine-side converter controllers of the high-speed FESS, which are designed for the applications of voltage and frequency support in low voltage grids. Description of the notations and values of parameters for the grid-side and machine-side converter controllers are provided in Table 5.3 and Table 5.4, respectively.

The grid-side converter controller has two major tasks:


**Grid voltage support:** Since the FESS in this work is connected to the low voltage distribution grid, it has to comply with the latest grid connection codes for this voltage level. Since November 2018, the new version of the grid code VDE-AR-N 4105 [182] is applicable for the German grids, which indicates that any ESS larger than 4.6 kVA, should follow the Q(U) characteristics, depicted in Figure 4.5(b). The Q(U) characteristics requires the ESS to start absorbing or injecting reactive power when the voltage passes the ±3 % of the nominal value and reaches its maximum reactive power capability when the voltage reaches ±7 %. The gradient of the Q(U) characteristics determines the voltage droop coefficient Dq, which indicates how much reactive power should be injected or absorbed with respect to the voltage changes. The voltage deadband udb is set to ±3 %, which helps avoiding continuous operation of the FESS during minor voltage variations. The Q(U) characteristic replaces the grid code requirement from 2011, which indicated that reactive power support is only required if the ESS injects power larger than 50 % of its rated value to the grid.

In this work, the instantaneous grid voltage (<sup>g</sup> ) is calculated in the dq-rotating reference frame and is used for the calculation of the reactive power reference, which enables an immediate detection of any changes in the voltage. This is calculated using

$$u\_{\rm g} = \sqrt{u\_{\rm dg}^2 + u\_{\rm qg}^2} \tag{5.20}$$

The reactive power reference for grid voltage support is generated according to the value of <sup>g</sup> , the voltage deadband udb, and the voltage droop coefficient Dq. The reference value for the reactive power can also sent by an operator using cref (see Figure 5.8).

**DC-link voltage control:** A conventional PI controller with saturation is used for controlling the DC-link voltage. The output of the PI controller determines the active power reference for the grid-side converter controller (gref), which is limited to the maximum power of the FESS (max).

**Determining the current reference values:** The upper level controllers or commands by the operator determine the active and reactive power reference values (gref and gref). The instantaneous power theory [237] is used for calculating the d- and q-axis current reference values (qgref and dgref) from the active and reactive power set points. Similar to the PMSM, the active and reactive power of the grid-side converter can be calculated from the voltages and currents in the dq-rotating reference frame using

$$P\_{\rm g} = \frac{3}{2} \{ \boldsymbol{u}\_{\rm dg} \boldsymbol{i}\_{\rm dg} + \boldsymbol{u}\_{\rm qg} \boldsymbol{i}\_{\rm qg} \},\tag{5.21}$$

and

$$Q\_{\rm g} = \frac{3}{2} \left( -u\_{\rm dg}i\_{\rm qg} + u\_{\rm qg}i\_{\rm dg} \right). \tag{5.22}$$

The dq-rotating reference frame of the grid-side converter controller is synchronized by the phase angle of the grid voltage (PLL), which is measured by a conventional Synchronous Reference Frame Phase-locked Loop (SFR-PLL). It is assumed that the dqrotating reference frame is chosen in a way, that the grid voltage is fully aligned with the d-axis, resulting in a q-axis component of zero for the voltage, i.e., qg = 0. The phaselocked loop controls the rotating reference frame angle, in order to maintain this condition with the changes in the system operating points (see subsection 5.2.5). By doing so, according to Eq. (5.21) and (5.22), the active and reactive power can then be independently controlled by controlling the values of dg and qg, respectively. Therefore, the d- and q-axis reference currents are calculated using

$$d\_{\rm dgref} = \frac{2P\_{\rm gref}}{3u\_{\rm dg}},\tag{5.23}$$

and

$$d\_{\rm qgref} = \frac{-2Q\_{\rm gref}}{3u\_{\rm dg}}.\tag{5.24}$$

The reference values for the d- and q-axis currents are sent to the inner current controller. The inner current controllers are conventional PI controllers. As seen in Figure 5.8, the outputs of the PI current controllers are added with decoupling terms in order to enable an independent control of the d- and q-axis currents. A feedforward term has also been added to reduce the high transient currents during the converter start-up, which also decouples dynamics of the converter system from those of the grid, and improves its disturbance rejection capability, as suggested in [243]. At the end, the grid-side converter controller generates the modulation index dg and qg in the dq-rotating reference frame, which are transformed back to the stationary frame using the inverse Park transform and the grid phase angle measured by the phased-locked loop.


Table 5.3: Variables definition and parameters of the grid-side converter (GSC) controller.

(a) Block diagram of the Grid-side Converter (GSC) controller. The upper part of the controller is responsible for supporting the grid voltage, while the lower part controls the DC-link voltage. (b) Block diagram of the Machine-Side Converter (MSC) controller. The upper part is responsible for the providing frequency support to the grid. The Maximum Torque per Ampere (MTPA) control is implemented for controlling the PMSM.

#### **5.2.4 Machine-side Converter Controller**

The machine-side converter controller controls the active power flow between the PMSM and the DC-link capacitor through controlling the currents of the machine. The active power reference can either be set by the grid frequency support block or by the operator through cref. The losses and auxiliary power required for running the entire system, represented by L&A in Figure 5.8(b), should also be considered for providing a certain power at the grid connection point. This is due to the fact, that this particular FESS has no additional auxiliary power connection, and all auxiliary components are also powered from its main terminal.

The active power of the flywheel is transferred to the grid by the grid-side converter controller, as it attempts to maintain the DC-link voltage. A description of the notations and values of parameters of the machine-side converter is provided in Table 5.4.

**Grid frequency support:** Similar to the voltage support requirements, the FESS has to follow the frequency support requirement of the latest German grid code for low voltage grids, the VDE-AR-N 4105:2018-11 [182]. This grid code indicates that the FESS has to stay connected to the grid in the frequency range of 47.5 Hz to 51.5 Hz, which is also known as the frequency ride-through. In this range, the FESS has to follow the P(f) characteristic shown in Figure 4.5(a). According to the P(f) characteristic the frequency deadband (fdb) for systems connected to the low voltage grids is 200 mHz. Also, the grid code requires different frequency droop coefficients for over- and under-frequency events larger than the deadband. Hence, the droop for the over-frequency (Dup) and underfrequency events (Ddn) are separated in Figure 5.8. In case of an over-frequency disturbance in the range of 50.2 to 51.5 Hz, the FESS has to participate in the frequency regulation by absorbing power with the rate of 40 % per Hz, after deducting the frequency deadband. In case of under-frequency events, the FESS has to inject power to the grid with the rate of 100 % per Hz until reaching its maximum power. Therefore, in Figure 5.8(b),

$$P\_{\rm mref} = \begin{cases} 20 \text{P\_n} \frac{50.2 - f\_{\rm g}}{50 \text{ Hz}} & \text{for } 50.2 \text{ Hz} \le f\_{\rm g} \le 51.5 \text{ Hz} \\ 50 \text{P\_n} \frac{49.8 - f\_{\rm g}}{50 \text{ Hz}} & \text{for } 47.5 \text{ Hz} \le f\_{\rm g} \le 49.8 \text{ Hz} \end{cases}' (5.25)$$

where P<sup>n</sup> is the nominal power of the FESS and <sup>g</sup> is the grid frequency.

**Speed-dependent ramp-rate limiter:** In practice, a FESS can only deliver the nominal power above a particular minimum rotational speed. This is due to the fact, that at low rotational speeds, providing the nominal power requires high amounts of torque, which can be above the maximum permissible torque of the PMSM. Therefore, it is important to consider the torque limiter in modeling a FESS, as shown in Figure 5.8(b). For this particular FESS, the torque limiter limits the power ramp rate to approximately 410 W/s, until reaching the speed of 23,200 rpm.

**Determining the current reference values:** There are several methods for controlling a PMSM [231], [235], [244], [247]. A categorization of some of the most well-known methods is shown in Figure 5.9. Selecting the most appropriate control method depends on the structure and type of the PMSM, and its main application. Conventional Field-Oriented Control (FOC) methods are generally more preferred due to a simpler design.

Figure 5.9: Different control strategies for controlling a PMSM (figure adopted from [247]).

For the surface mounted PMSM used in the high-speed FESS, the Maximum Toque per Ampere (MTPA) control has been used. This method is suitable for applications in which fast reaction times are required. With the use of the maximum torque per ampere control, the copper losses of the machine is also minimized [248]. As shown in Eq. (5.8), in a surface mounted PMSM, the electrical torque is a function of the q-axis component of the stator current only. Therefore, for a given stator current, the maximum torque is achieved, when the current is fully aligned with the q-axis. Therefore, to implement the maximum torque per ampere control, the d-axis component of the current should be forced to zero, i.e., dmref = 0. Figure 5.10(a) shows the vector diagram of the stator current without the maximum torque per ampere control, while Figure 5.10(b) depicts the case when the maximum torque per ampere is applied. This method is also known as the constant torque angle control, since the torque angle, in Figure 5.10, is always kept at 90 degrees with respect to the flux linkage of the permanent magnets (PM). The reference value for the q-axis current (qmref) is calculated from Eq. (5.8) from the reference value for the electrical torque (eref) using

$$d\_{\rm qmref} = \frac{2\tau\_{\rm eref}}{3n\_{\rm p}\psi\_{\rm PM}}.\tag{5.26}$$

The torque reference itself is calculated from the active power reference of the machineside converter controller (mref) and the instantaneous speed of the FESS (m). Conventional PI controllers with decoupling terms have been used for the inner current controllers, similar to the grid-side converter.

Table 5.4: Variables definition and parameters of the machine-side converter (MSC) controller.


It is important to note that in high-speed machines, such as the PMSM in a high-speed FESS, the speed and position of the rotor are often not directly measured, and are estimated from the voltage and current values, with the methods described in [249], [250]. However, since these estimation methods often show a good accuracy at high speeds [249], in which the FESS is often operated, the estimation procedure is neglected in this work, and rotor speed is directly calculated from the machine model. This is also necessary to limit the model computation time for real-time simulations.

### **5.2.5 Modeling of Synchronous Reference Frame Phase-Locked Loop**

For some applications, including frequency support in low voltage and low inertia grids (see Chapter 7), the dynamics of the phased-locked loop should also be modelled [110]. In this chapter, a Synchronous Reference Frame Phase-locked Loop (SRF-PLL) with a low-pass filter, as shown in Figure 5.11, is used for tracking the frequency and phase of the grid voltage. The estimated phase angle by the phase-locked loop (PLL) is used in the grid-side converter controller to synchronize the dq-rotating reference frame with the grid voltage. This is also necessary for connecting the FESS to the grid with minimum transient currents.

As seen in Figure 5.11, the grid voltage, after being converted to its per unit value, is multiplied by the cosine of the phase estimated by the phased-locked loop. By assuming that the grid voltage in per unit is sin( + ), this results in

$$\sin(\omega t + \varphi)\cos(\theta\_{\rm PL}) = \frac{1}{2} \left( \underbrace{\sin(\omega t + \varphi + \theta\_{\rm PL})}\_{\text{Periodic component at } 2\omega} + \underbrace{\sin(\omega t + \varphi - \theta\_{\rm PL})}\_{\text{DC component}} \right). \tag{5.27}$$

If the phase estimated by the phase-locked loop is approximately close to the phase of the input signal, i.e., PLL ≈ + , then Eq. (5.27) consists of a term with double the fundamental frequency and a DC component close to zero. This is due to the fact, that sin( + − PLL) in case of small phase angles is approximated by + − PLL. The moving average block in Figure 5.11 eliminates the periodic component by averaging at the fundamental frequency. Therefore, only the DC component remains, which reflects the error in estimating the phase of the input signal by the phase-locked loop. A PI controller forces the steady-state error to zero. The output of the PI controller is the electrical angular velocity. A simple numerical integration generates the phase output of the phase-locked loop. The angular velocity is then converted to Hertz and passed through a low-pass filter and rate limiter to reduce the effect of harmonics and asymmetries of the input voltage on the estimated frequency [251]. The PI controller parameters used for the phase-locked loop in this work are tuned to have a quick response to step changes in the frequency, and are given in Table 5.5.

Figure 5.11: Block diagram of a synchronous reference frame phase-locked loop with a low-pass filter.

Table 5.5: Variables definition and parameters of the synchronous reference frame phase-locked loop (SRF-PLL).


# **5.2.6 Modeling of System Losses and Auxiliary Power Requirements**

For obtaining an accurate model of a FESS, it is also necessary to model the losses within the system. The losses in a FESS can be categorized in three main groups, which are:


the copper losses in the PMSM, the two voltage source converters, filters, and other conductors within the FESS.

▪ **Constant losses:** This includes any losses within the system components that are constant throughout the whole operation range of the FESS.

As discussed, the high-speed FESS in this work has only one grid connection point. Therefore, the power required for all the axillary components including the cooling water and vacuum pumps, active magnetic bearings, control and measurement systems, fans and container ventilation system, fire and safety systems, and all other components are also provided by the flywheel. The auxiliary power is constant for some components, such as the power supply for the controllers, while for others, such as for the cooling systems, can vary depending on the system operating conditions, including the speed and the current of the FESS.

Due to the complex structure of the high-speed FESS with numerous components, the losses and auxiliary power are aggregated all together and estimated directly from measurement, rather than using an analytical approach to calculate each part individually. The three types of losses in the FESS and the axillary power required for running the system are aggregated and modeled in one variable, L&A, which is calculated using

$$P\_{\rm L\&A} = c\_1 + c\_2 l\_s^2 + c\_3 \mu\_{\rm m}^2. \tag{5.28}$$

In Eq. (5.28), <sup>1</sup> is the constant component for constant losses and auxiliary power requirements, <sup>2</sup> is the coefficient for current-dependent losses, while <sup>s</sup> is the current at the FESS terminals, and <sup>3</sup> is the coefficient for speed-dependent losses. The parameters 1 , <sup>2</sup> , and <sup>3</sup> are estimated using real measurements from the FESS. The experimental setup and measurements are described later in Chapter 6. A simple curve fitting approach based on least-square error method is used to estimate the parameters of Eq. (5.28), and their values are given in Table 5.6. These values are estimated based on the difference between the power that the flywheel delivers and the power at the FESS terminals. It can be seen from Table 5.6 that the constant component of Eq. (5.28), which mainly consists of the auxiliary power required to run the system, is significant, considering a 60 kW FESS. However, it should be mentioned that some of the auxiliary components are meant for multiple flywheels, and the required auxiliary power for running multiple flywheels does not significantly increase with the number of flywheels.


Table 5.6: Parameters of the total system losses and auxiliary power requirements.

### **5.3 Real-time Simulation of the High-speed FESS**

In this section, the model of the high-speed FESS, described in section 5.2 using the model of each of its component, is simulated in real-time, and the results are compared with offline simulation results. For the real-time simulations, the model has been implemented in Opal-RT's RT-LAB real-time simulation environment. The real-time simulations are carried out on the Opal-RT's 5600 digital real-time simulator. RT-LAB generates a precompiled C code from the model on a host PC and loads this code on the real-time target. For offline non-real-time simulations, the ordinary differential equations of each component of the FESS, given in section 5.2, are discretized and solved using numerical methods in a script written in MATLAB programming language. The same simulation time step of 10 µs has been used in both cases. The simplified power system shown in Figure 5.12 is used to demonstrate the performance of the FESS. The equivalent grid and the transformer impedance are based on the values from the CIGRE European low voltage benchmark grid, in which a short-circuit capacity of 100 MVA and an X/R ratio of 1 has been considered for the medium voltage grid [123]. The load is a simple RL load. The model comparison is made using steps in the grid frequency and voltage.

Figure 5.12: The single-line diagram of the simulated system.

#### **5.3.1 Frequency Step Response**

To show the fast dynamic behavior and full response of the FESS, a step of 1.2 Hz in the grid frequency is simulated, from 50 Hz to 48.8 Hz to fully observe the flywheel dynamics. Since the FESS is connected to a stiff grid, it cannot influence the frequency. However, it should provide active power to the grid according to its droop characteristic. In this case, the FESS should inject 100 % of its nominal power to the grid, according to the latest German grid code [182]. The high-speed FESS is initially charged to reach 36,000 rpm, with no power being exchanged with the grid prior to the simulated event.

The active power of the FESS model in real-time and offline simulation is shown in Figure 5.13, while other variables of the FESS are shown in Figure 5.14. The step in the frequency reference is applied at 0.04 s after the start of the simulation. As seen, the FESS quickly reaches its nominal power of 60 kW, and its rotational speed, shown in Figure 5.14(b), and therefore, its energy content begins to fall. Since only a short time span of 0.16 s is shown here to show the fast system dynamics, the drop in the flywheel speed is not significant. The DC-link voltage, shown in Figure 5.14(c), is quickly restored to 720 V by transferring the power injected from the PMSM to the grid by the grid-side converter. The electrical torque of the machine and the q-axis current required for providing this torque are depicted in Figure 5.14(d) and (f), respectively. As seen, a good match between the real-time simulation results obtained from the RT-LAB simulation software and the offline non-real-time results of the MATLAB script can be observed in all FESS variables. Two other steps in the grid frequency with amplitudes of 49.3 Hz and 49.75 Hz have also been applied for testing the response of the FESS model. The active power of the FESS from the real-time and non-real-time simulations for these steps in the frequency are also shown in Figure 5.13. It can be seen that in all cases the real-time simulation results in RT-LAB match the offline results for the MATLAB script.

Figure 5.13: The active power injected by the FESS during different steps in the grid frequency

Figure 5.14: Real-time and non-real-time simulation results during a step in the grid frequency by 1.2 Hz.

#### **5.3.2 Voltage Sag Response**

In this scenario, a voltage sag on the medium voltage grid is simulated. Such voltage drops are common and could be caused by faults in the high voltage grids or parallel feeders [254]. The depth of the voltage sag is assumed be 7 % of the nominal values. During such a disturbance, the FESS should inject its maximum reactive power to the grid. The real-time and offline simulation results for this scenario are shown in Figure 5.15. The grid-side converter of the FESS injects approximately 50 kVar of reactive power, which the converter's maximum reactive power, to support the load during the voltage disturbance, as seen in Figure 5.15(a). This is done by injecting current in the qaxis direction, which is shown in Figure 5.15(e), as the reactive power is controlled by qaxis current of the grid-side converter. Since only reactive power is provided by the FESS, there is almost no change in the FESS speed, shown in Figure 5.15(b), the DC-link voltage, shown in Figure 5.15(c), and other variables of the machine and machine-side converter. Some minor transients are observed only at the beginning of the voltage sag. The speed of the FESS is slightly decreasing due to the losses within the FESS and the auxiliary power requirement. More importantly, the real-time simulation results match the non-real-time offline simulation in both steady-state and transient conditions for all variables of the FESS.

Figure 5.15: Real-time and non-real-time simulation results during a 7% voltage sag.

# **6 PHIL Testing and Model Validation of a High-Speed FESS for Frequency Support**

This chapter presents the Power Hardware-in-the-Loop (PHIL) testing results of a 60 kW high-speed FESS for the application of frequency support. Firstly, the characteristics of the FESS, as the hardware under test, are given, followed by a description of the PHIL testing facility built and used for verifying the performance of the FESS. Next, several frequency deviation scenarios are simulated in real-time in order to evaluate the response of the FESS during such disturbances. The PHIL testing results are also used to validate the model of the high-speed FESS, described in Chapter 5. The same inputs which have been sent to the real FESS during the PHIL tests, are also given to the model of the FESS, running in real-time, and its response is compared with the response of the real FESS.

As discussed in Chapter 3, no PHIL testing of a FESS has been reported previously for verifying its response during frequency disturbances. There are only some limited experiments carried out for the application of supporting pulsed power loads in DC grids of shipboard power systems [24], [25]. Regarding the model validation of a FESS, an experimental validation of a FESS model has been reported in [255]. However, the tested system is a laboratory-scale 5.5 kW low-speed FESS with a maximum rotational speed of 3000 rpm, which is not a practical system with adequate capacity and power that can be used for ancillary services in the grid. In addition, no frequency control loop was implemented and tested. Therefore, this chapter presents the first PHIL testing of a highspeed FESS for frequency support, and the first model validation of a high-speed FESS.

## **6.1 Description of the Hardware under Test**

The hardware under test is the 60 kW high-speed FESS, which has been modeled in Chapter 5. The inside view of the container of the FESS with all the main and auxiliary components is shown in Figure 5.2. This FESS has a maximum rotational speed of 45,000 rpm, which corresponds to a maximum capacity of 3.6 kWh, taking into account the flywheel's inertia. Therefore, this FESS can provide a discharge time of several minutes, depending on the discharge power. The high rotational speed of the FESS is achieved with the use of a carbon fiber reinforced plastic rotor, which is rotated using active magnetic bearings in a vacuum enclosure. The vacuum pressure of the flywheel

containment is below 0.03 mbar during normal operation of the system. The flywheel steel containment is mechanically locked to a large underground cement foundation for safety purposes.

# **6.2 Description of the PHIL Setup**

The PHIL testing of the FESS is conducted at the KIT's 1 MVA PHIL testing facility, as part of the EnergyLab 2.0, which is described in section 3.3. A schematic diagram of the setup tailored for the PHIL testing of the FESS is shown in Figure 6.1.

Figure 6.1: The schematic diagram of the setup used for the PHIL testing of the 60 kW high-speed FESS.

An Opal-RT 5700 real-time simulator is used for simulating the grid models and scenarios in real-time, controlling the power amplifier, collecting and recording measurement data, and monitoring and controlling the FESS operation during the experiments. The grid simulations, data collection, and the data conversion required for the Small Form-factor Pluggable (SFP) connection to the amplifier (see subsection 3.3) is carried out on one CPU core. A separate core of the simulator is dedicated for an Ethernet connection between the Programmable Logic Controller (PLC) of the FESS and the realtime simulator using the Modbus TCP/IP protocol. This is used to read the internal variables from the FESS controller and to send control commands to it. A list of the most important commands and measurement values communicated through the Ethernet connection with the FESS controller is given in Table 6.1.

Table 6.1: List of the most important commands and measurement values communicated through the Ethernet connection using the Modbus TCP/IP protocol between the FESS controller and the realtime simulator.


A single 200 kVA GAMP6 switched-mode power amplifier from Egston Power company is used for energizing the FESS using a four-wire three-phase cable. The power amplifier receives the voltage set points from the real-time simulator using a dedicated SFP connection. The voltages and currents of the FESS are measured at the FESS terminals using voltage and current transducers, and are sent to two separate Opal-RT OP4520 Kintex7 FPGA & I/O expansion units, due to the long physical distance between the hardware under test and the real-time simulator. Current and voltage measurements have been recalibrated before the experiments using a precise power supply.

Before and after the experiments, the vacuum, cooling water pumps, and several other auxiliary components of the FESS have to run continuously. Therefore, a manual switch was installed at the FESS terminals, which can switch the supply of the FESS from the power amplifiers to a conventional power supply, when no experiments are conducted.

The voltage-type ideal transformer method [256], which has been described in detail in subsection 3.2.2.1, is used for interfacing the hardware under test and the real-time simulator. For the current feedbacks a first-order low-pass filter with a cut-off frequency of 100 kHz is applied. A stable PHIL setup was observed using the ideal transformer method during all experiments and scenarios. Therefore, there has been no need to apply more complex PHIL interfacing algorithms or reduce the cut-off frequency of the filter.

# **6.3 PHIL Testing Scenarios and Model Validation of the FESS**

The high-speed FESS is tested in the following scenarios:


The voltages, commands, and other input signals sent to the real FESS during the PHIL tests in each scenario are also given to the FESS model, running in real-time, for model validation. In the following, the PHIL testing results together with the response of the model of the FESS are provided for each scenario, and the response of the FESS and its model are compared with each other.

It should be noted that in all cases the FESS is tested in a closed loop with the current feedbacks being injected into the real-time simulation.

### **6.3.1 Scenario A: Charging at the Rated Power**

This scenario includes charging the FESS from being completely discharged at a standstill condition to reaching its maximum rotational speed and the nominal capacity. The aim of this test is to derive the charging profile of the FESS, in particular at low rotational speeds, where the torque limitations of the machine prohibit the FESS from absorbing the nominal active power.

A simplified power system consisting of a voltage source and equivalent impedance with the values from the CIGRE European low voltage grid benchmark [123] is simulated in real-time on the Opal-RT's OP5700 real-time simulator. The voltages after the equivalent impedance are sent to the power amplifier, which feeds the FESS. Next, a step in the active power reference of the FESS from 0 to -60 kW is applied from the OP5700 realtime simulator using the Ethernet connection with the FESS controller and the Modbus protocol. The FESS should begin charging itself. When the FESS is fully charged and reached a state of charge of 100 %, which corresponds to a rotational speed of 45,000 rpm for this particular FESS, it should automatically stop absorbing power.

The results obtained from the PHIL testing of the FESS and the real-time simulation of its model in this scenario are depicted in Figures 6.2 and 6.3. The active power reference of -60 kW is applied at 5.7 s after the start of the test. As expected, the FESS can only absorb the nominal power above a certain rotational speed, due to the machine's torque limitations. This can clearly be seen in Figure 6.2, in which the ramp rate of the active power is limited, until reaching a state of charge of 24.6 %, which corresponds to a rotational speed of approximately 23,200 rpm. After exceeding this point, the FESS can absorb the rated active power until being fully charged. At the end, when the FESS is fully charged, the active power of the FESS drops to the minimum power required to cover the auxiliary power and losses in the whole system, as seen in Figure 6.2. This minimum power is required to maintain the state of charge at the maximum value.

It can be seen from Figures 6.2 and 6.3 that the real-time simulation model can reflect the charging profile of the real FESS with an acceptable accuracy. The time required to fully charge the FESS is approximately the same in both simulation and experimental results and the variations in the state of charge of the FESS during the charging are very close, with a maximum difference of just below 0.8 %. This small error could be due to the variations in the auxiliary power of the FESS during the experiment.

Figure 6.2: The measured and simulated active power of the FESS during the PHIL testing of scenario A, charging the FESS at the rated power.

Figure 6.3: The measured and simulated state of charge of the FESS during the PHIL testing of scenario A, charging the FESS at the rated power.

## **6.3.2 Scenario B: Grid Code Compliance Verification Test for Frequency Support**

The aim of this test is to verify that the high-speed FESS follows the P(f) characteristics during frequency deviations, according to the latest German grid code for connections to the low voltage distribution grids, the VDE-AR-N 4105:2018-11 [182]. The P(f) characteristic, shown in Figure 4.5(a), indicates the frequency deadband and different droop coefficients for over- and under-frequency incidents. According to this grid code, for frequency deviations above 200 mHz, the FESS should react within 2 s and reach the active power set point given by the P(f) characteristic with a rise time of under 2 s and a settling time of under 20 s. The acceptance tolerance for the difference between the active power set point and the measured active power is ±10 % [257].

In order to verify the grid code compliance of the FESS, a simplified power system similar to the previous scenario is simulated in real-time on the OP5700 real-time simulator. After charging the FESS, the frequency of the simulated grid is altered according to Figure 6.4. This test is adopted from the VDE V 0124-100:2020-06 [257], which describes the testing protocols for the generation units and storage systems that are planned to be connected to low voltage distribution grids. Conducting these tests is now a mandatory requirement for any ESS or generation unit, prior to a grid connection [257]. In comparison to the official testing requirements of the grid code, the tests have been shortened, taking into account the limited energy content of the FESS. It can be seen from Figure 6.4, that this test includes steps of different amplitudes in the emulated grid frequency. Each step in the frequency is kept for the period of 30 s, after which the frequency is set back to the nominal value. Each over-frequency step is followed by a subsequent under-frequency step with the same amplitude.

Figure 6.4: The variations in the grid frequency for grid code compliance validation of the high-speed FESS, adopted from the grid code VDE V 0124-100:2020-06 [257].

The tests are conducted with the FESS operating at the power factor of unity, as requested by the testing protocol [257]. The test is carried out after charging the FESS to reach a state of charge of 75 %, which corresponds to approximately 39,000 rpm for this FESS. At this state of charge, there are no speed-dependent ramp rate limitations due to the machine's maximum torque, and the FESS will provide the highest possible ramp rate, taking into account other constraints of the system, such as the measurement and communication delays.

The PHIL testing results of the 60 kW FESS and the real-time simulation results for this scenario are shown in Figures 6.5-6.7. The measurements from the FESS are obtained from the FESS internal controller, which shows some discontinuity in the results, due to the use of the Modbus protocol. It is clear from Figure 6.5 that the FESS controller can easily track the changes in the grid frequency with a high accuracy during frequency deviations. More importantly, it is shown in Figure 6.6, that the FESS quickly reacts to frequency changes larger than the frequency deadband, while ignoring smaller changes in

the frequency. In the case of over-frequency steps, the FESS absorbs power from the grid with the rate of 40 % per Hz, after deducting the frequency deadband. For instance, it can be seen that the FESS absorbs 24 kW of power, when the emulated grid frequency is set to 51.2 Hz. During an under-frequency event, the FESS injects power to the grid with the rate of 100 % per Hz, after deducting the frequency deadband. As an example, it can be seen that when the simulated grid frequency is set to 49.3 Hz, the FESS injects 30 kW to the grid, which is 50 % of its nominal power. Therefore, it can be concluded the FESS accurately follows the P(f) characteristic required by the German grid code for an ESS connected at the low voltage level.

The simulation model of the FESS also follows the P(f) characteristics of the grid code, and a good match between the simulated power and the measured power of the FESS is observed in Figure 6.6. Therefore, the model of the FESS can accurately generate the same output of the real hardware. This is also shown in Figure 6.7, where the state of charge of the FESS is depicted. Between the state of charge of the real system and the simulation model, a maximum mismatch of only 0.35 % is observed in this scenario.

Furthermore, the fact that the FESS remains connected throughout these extreme frequency deviations and is not automatically disconnected from the grid by its protection systems demonstrates the frequency ride-through capability of the FESS according to the requirements of VDE-AR-N 4105:2018-11 [182].

Figure 6.5: The measured and simulated grid frequency during the PHIL testing of scenario B, grid code compliance verification test.

Figure 6.6: The measured and simulated active power of the FESS during the PHIL testing of scenario B, grid code compliance verification test.

Figure 6.7: The measured and simulated state of charge of the FESS during the PHIL testing of scenario B, grid code compliance verification test.

In order to see how fast the real FESS responds to the grid frequency changes, one of the emulated steps in the grid frequency is shown more closely in Figure 6.8. Here, the active power is calculated directly from the voltage and current measurements at the FESS terminals, since the measurements from the FESS inertial controller have a high latency. It can be seen from Figure 6.8 that the FESS reaches the required active power in 58.5 ms, and settles in 88.6 ms. This is a very fast response time from the FESS. It is important to note that this time includes the time required for the phase-locked loop of the FESS to measure the grid frequency, the delay in signal transmission, and more importantly, the delay caused by the filter applied on the active power measurement. Therefore, the actual response time of the FESS can be even faster. This is by far faster than the requirements of the grid code in Germany, which requires an ESS to react to the frequency incident in below 2 s, to reach the set point in another 1 s, and to settle in below 20 s [257]. This delay is an important parameter, which is commonly neglected in studies regarding fast frequency support of converter-interfaced systems in power systems [12].

Figure 6.8: The response time of the high-speed FESS to a step in the grid frequency during the PHIL testing.

## **6.3.3 Scenario C: Response to the UK's Frequency incident of August 9, 2019**

In this scenario, the recorded frequency variations during the major grid disturbance of August 9, 2019, in the power grid of the UK is simulated in real-time in order to evaluate the response of the high-speed FESS during such frequency deviations.

On this day, the UK's power system experienced a severe frequency drop, which led to the loss of power for approximately 1.1 million customers, including hospitals, airports, and railways [258]. The sequence of events, which is shown in Figure 6.9, was triggered by a lightning strike on a transmission line. The fault caused by the lightning strike resulted in a generation loss of approximately 1131 MW, which included 150 MW of distributed energy resources connected to the low voltage distribution grids, and reduction in the power generated by a wind farm by 737 MW. This significant loss in generation led to a severe drop in the grid frequency, reaching 49.1 Hz. As the frequency began to recover, another 210 MW gas turbine tripped, followed by an additional 350 MW power loss due to the activation of the Rate of Change of Frequency (ROCOF) protection system of several synchronous generators. The frequency dropped again, reaching a frequency nadir of 48.8 Hz. At this point, an under-frequency load shedding was enabled, disconnecting approximately 5 % of the total load, which included roughly 1000 MW of the loads at the low voltage distribution grids. With the activation of the under-frequency load shedding scheme, the frequency was regulated back to 50 Hz. The changes in the grid frequency in the power system of the UK during the time span of these incidents, which lasted approximately around 5 minutes, are shown with the blue dotted line in Figure 6.10.

Figure 6.9: Sequence of event involved in UK's major frequency incident of 9 August, 2019 [258].

The original measured frequency data from the National Grid of UK has been used for testing the response of the FESS [259]. A simplified power system similar to the previous systems is simulated in real-time, where the frequency can follow the recorded frequency of the incident, when required. The FESS droop and deadband settings of the FESS are kept according to the P(f) characteristics of the German grid code VDE-AR-N 4105:2018-11 [182]. Prior to activating the frequency disturbance, the FESS is charged to reach a state of charge of approximately 95 %, which corresponds to the rotational speed of approximately 43,860 rpm.

The PHIL testing results of the 60 kW high-speed FESS and the real-time simulation results of its model during this scenario are shown in Figures 6.10-6.12. It can be observed from Figure 6.10, that similar to the previous scenario, the controller of the high-speed FESS easily tracks the changes in the grid frequency. Also, the FESS injects power to the grid, as soon as the frequency violates the frequency deadband of 200 mHz, and according to its droop setting, as seen in Figure 6.11. Since the maximum frequency deviation is approximately 1 Hz after the frequency deadband, the FESS will reach its maximum power of 60 kW at the minimum point of the frequency.

In addition, an excellent match between the active power of the real FESS and its model is observed in Figure 6.11, which indicates that both the real hardware and its model can successfully track the active power reference generated from the P(f) characteristics of the grid code. More importantly, the state of charge of the real hardware and its model during the time period of the disturbance are compared with each other in Figure 6.12. It can be seen, that there is a good match between the state of charge of the real FESS and the model in this scenario as well. There is a slight difference between the state of charge values of maximum 0.7 %, which could have been caused by the variability of the auxiliary power of the real FESS, such the power required for running the pumps and the cooling systems. It is also interesting to observe that the high-speed FESS has enough capacity to provide the required support throughout the whole frequency incident.

Figure 6.10: The measured and simulated grid frequency during the PHIL testing of scenario C, response to the frequency incident of August 9, 2019 in the UK.

Figure 6.11: The measured and simulated active power of the FESS during the PHIL testing of scenario C, response to the frequency incident of August 9, 2019 in the UK.

Figure 6.12: The measured and simulated state of charge of the FESS during the PHIL testing of scenario C, response to the frequency incident of August 9, 2019 in the UK.

# **6.3.4 Scenario D: Increasing the Frequency Support of a Low-voltage Distribution Grid**

The goal of this scenario is to show the increase in the contribution of a low voltage distribution grid for supporting the grid frequency by using a FESS during frequency incidents.

A real low voltage distribution network in southern Germany is modeled and simulated in real-time in a collaboration with the distribution system operator. The single-line diagram of the modeled grid is depicted in Figure 6.13. This low voltage grid includes 78 short underground cables, 35 load connection points, and 3 PV systems. One of the PV systems is a medium-scale 100 kW PV system, while others are small residential ones with ratings below 10 kW. The medium voltage grid has an equivalent short-circuit capacity of 42.2 MVA with an X/R ratio of 1.61. These values are calculated based on the model of the medium voltage grid provided by the system operator, using the network reduction tool and extended Ward method [260] in DIgSILENT PowerFactory. The transformer and line data are also provided by the system operator.

Figure 6.13: The single-line diagram of the real low voltage distribution network in southern Germany, used for PHIL testing of the FESS in scenario D.

The PV systems are modeled based on the generic two-stage PV model described in [261]. Conventional control loops including frequency control, active and reactive power control, DC-link voltage control, and maximum power point tracking using the incremental conductance method are also included in the PV models. For the PV cells, the parameters from of a commercial PV cell type, the Solarwatt Vision 60M style [262] are used.

Using conventional power system solvers such as state-space solvers, this network cannot be simulated in real-time with a simulation time step of below 200 µs, even without considering the PV systems. Therefore, it is not feasible to simulate accurately this grid model in real-time. This is due to the large number of short three-phase underground cables. Since the longest cable in this model is only about 300 m long, the conventional model decoupling method at long lines cannot be applied here, and one or several of the techniques listed in Table 3.1 have to be applied to be able to run this model in real-time.

Table 6.2 shows the required execution time for this grid model using different techniques listed in Table 3.1 to reduce the model execution time, which were calculated using the RT-LAB software. Firstly, the State-space Nodal (SSN) solver [109] is used to decouple the system equations into 12 subgroups, referred to as the nodal groups. This technique adds no artificial delays or impedance changes to the model, so it does not compromise the model accuracy. As shown in Figure 6.13, the secondary side of the transformer is selected as the main decoupling point (SSN node), as it generates a high number of groups from a single point. Using the state-space nodal solver, this grid model without the PV systems can be simulated in real-time with a simulation time step of 40 µs. Another method to reduce the model execution time is to neglect the line capacitances. This can be done for this particular grid as lines are mostly below 100 m long, are operated at the low voltage level, and the grid is assumed to be symmetrical [117]. This reduces the model execution time of the grid model without the PV systems to 20 µs, as it reduces the number of state variables by 8 for each line, and the grid under study has a relatively large number of lines. When compared to the model with the line capacitances, the grid model without the line capacitances shows on average a 0.06 % error in the steady-state currents and 0.14 % error in the steady-state voltages. By combining the use of state-space nodal solver and neglecting the line capacitances, the required model execution time for the grid without the PV systems decreases to only 8.72 µs. To add the PV models, the average converter models are used, as the study focuses on the frequency control [110]. By adding the PV models, the model execution time increases to 17.6 µs. Considering the added components to the model required for the communication to the power amplifier and with the two I/O expansion units, and also software safety features to shut down the system when extreme values are observed during PHIL testing, the simulation time step is set to 24 µs, in order to avoid any unexpected overruns. This relatively low simulation time step leads to a small loop delay. This improves the stability of the PHIL setup, as discussed in subsection 3.2.2.

The real-time simulation results of the grid model using state-space nodal solver are validated against offline results in MATLAB/Simulink, and no noticeable difference is observed. Moreover, the steady-state results are validated with the load flow results provided by the distribution system operator in PowerFactory, and the average difference between the two models was only 0.23 % in the steady-state voltages, and 0.3 % in the steady-state currents.


Table 6.2: Comparison of the model execution times of the German low voltage distribution grid using different techniques used for reducing the model execution time.

The real high-speed FESS has been integrated into the grid model at the secondary side of the main transformer, as depicted in Figure 6.13. A separate nodal group is dedicated for the connection of the FESS using the current values measured at the real FESS terminals.

A slow frequency ramp is simulated in the medium voltage grid with the frequency slowly increasing from 50 Hz to 51.2 Hz with a rate of 0.25 Hz/min. The real FESS together with the simulated PV systems react to the frequency ramp according to the German grid code for low voltage grids. The grid code indicates that both PV generation and storage units have to participate in the over-frequency incidents and reduce their output with the rate of 40 % per Hz, when the frequency passes the 200 mHz threshold.

The FESS is charged to approximately 71 % before conducting this experiment. The results of the PHIL testing of the FESS together with the results of the real-time simulation of the grid in this scenario are shown in Figures 6.14 and 6.15. At 10 s after the beginning of the test, the PV systems are switched on with an enabled maximum power point tracking at the irradiance density of 1 kW/m<sup>2</sup> and an ambient temperature of 25 ֯C. Therefore, they quickly reach their maximum power output, as seen in Figure 6.14(a). This reduces the total power that low voltage grid draws from the medium voltage grid from 140 kW to approximately only 25 kW. Such low powers are common around noon in low voltage grids with a high share of PV generation, as previously seen in Chapter 4. The frequency ramp is enabled at 27.4 s. As soon as the frequency passes the frequency deadband, the simulated PV systems starts to curtail their generation during the over-frequency event, as seen in Figure 6.14(a), and contribute to decreasing the total generation in the system. At the same time, the real FESS starts to draw active power from the grid according to its droop coefficient and begins to charge itself, as shown in Figure 6.14(d). This increases the consumption in the low voltage grid. When the frequency reaches 51.2 Hz, the total active power of the low voltage grid has increased from 25 kW to approximately 94 kW (seen Figure 6.14(c)). Approximately one third of the increase in power is the share of the real FESS, while the rest is due to the curtailment of the generation in the simulated PV systems. Therefore, the active components of the low voltage grid contributed to the regulation of the frequency of the whole system.

Figure 6.14: The PHIL testing results of scenario D, increasing the frequency support of low-voltage grids. (a) Power of the PV systems, (b) the measured and simulated grid frequency, (c) The total power drawn by the low voltage (LV) grid, and (d) the measured and simulated power of the FESS.

Figure 6.15: The measured and simulated state of charge of the FESS during the PHIL testing of scenario D, increasing the frequency support of low-voltage grids.

In addition, similar to the previous scenarios, a good match between the active power of the real FESS and the simulated one, running in real-time, is observed, as seen in see Figure 6.14(d). Also, it is clear from Figure 6.15, that the model can reflect the variations in the state of charge of the FESS with a high accuracy in this scenario as well.

### **6.4 Summary**

This chapter presented the PHIL testing results of a 60 kW/3.6 kWh high-speed FESS in several frequency deviations scenarios. It was demonstrated the high-speed FESS can accurately follow the P(f) characteristics required by German grid code for low voltage grids during under- and over-frequency disturbances. Furthermore, it was shown that the FESS can reach the required power in under 60 ms during frequency disturbances and support the grid by quickly charging and discharging itself. Moreover, the nonlinear charging profile of the FESS from a standstill condition until reaching the maximum rotational speed was demonstrated, showing the limitations of a FESS at low rotational speeds.

In addition, it was shown that the model of the FESS, presented in Chapter 5, can approximately generate the same power and state of charge of the real FESS in real-time, during all the investigated scenarios.

# **7 Adaptive Inertia Emulation using High-speed Flywheel Energy Storage Systems**

As discussed in Chapter 1, the continuous growth in the share of converter-interfaced renewable generation has led to a steady decrease in the total inertia in power systems. The reduction in the system inertia can lead to an increase in the Rate of Change of Frequency (ROCOF) and the severity of the frequency deviations during imbalances in generation and demand. Inertia emulation techniques using converter-interfaced systems can help to decrease the effects of low physical inertia in the system by providing synthetic inertia support during frequency disturbances.

In this chapter, a new adaptive inertia emulation controller for high-speed FESS is introduced. The proposed controller aims to reduce simultaneously the ROCOF and the frequency nadir by adaptively changing the inertia and damping coefficients of the controller. It uses a combination of the bang-bang control approaches and self-adaptive ones to benefit from the advantages of both control designs. The performance of the proposed adaptive controller is validated initially using offline numerical simulations in a low voltage microgrid. Next, the proposed controller is implemented on the 60 kW highspeed FESS, described in chapter 5 and 6, using the concept of rapid control prototyping. The performance of the FESS using the proposed inertia emulation controller is validated using PHIL testing of the FESS in the low voltage microgrid in two different scenarios. The simulation and PHIL testing results confirm that the use of the FESS with the proposed inertia emulation controller can effectively reduce the ROCOF and the frequency deviation, and the suggested adaptive control design outperforms several adaptive methods previously reported in literature.

# **7.1 Introduction**

In this section, the need for inertia emulation using converter-interfaced systems is discussed briefly. The impact of the cumulative system inertia on the frequency dynamics of a power system during a power imbalance is demonstrated using a simplified singlemachine power system. A comparison among the possible active power sources required for emulating the inertia response using converter-interfaced systems is also presented, aiming to justify the advantages of using a FESS for this application. Lastly, a literature review on different adaptive control methods for providing inertia support is provided.

## **7.1.1 The Need for Inertia Support from Converterinterfaced Systems**

Inertia of a power system resists fast variations of its frequency during disturbances. Analyzing the previous major frequency disturbances in Europe reveals that a stable system operation can be guaranteed only up to a ROCOF of 1 Hz/s, which corresponds to a system imbalance ratio of 20 %. However, imbalance ratios of more than 40 % are expected for the near future in Europe, which leads to a ROCOF of more than 2 Hz/s [5]. In poorly interconnected grids, such as in Australia, a ROCOF as high as 6 Hz/s has already been experienced in 2016, which led to a widespread system blackout, due to the activation of under-frequency and ROCOF protection relays of the synchronous generators of conventional power plants [263].

Figure 7.1 illustrates a typical frequency behavior of a power system during a major loss in generation. The initial slope of the drop in the frequency, or the ROCOF, is limited only by the cumulative system inertia. The system inertia plays an important role in maintaining the grid stability [264]. The inertia also has an impact on the minimum frequency reached during the under-frequency event, which is called the frequency nadir [265]. A low-inertia power system can experience a high ROCOF and an extreme frequency nadir during disturbances.

Figure 7.1: Typical behavior of the frequency in a power system after a major large loss in generation (adopted from [266]).

The impact of the system inertia on the frequency and the ROCOF can be seen more clearly if we simplify the frequency control loops of a power system to the block diagram shown in Figure 7.2, as proposed in [267]. In Figure 7.2, a multi-machine system is represented by an equivalent single machine, where the inertia constant H<sup>g</sup> is the equivalent inertia constant of all generating units, and the damping coefficient D is the equivalent damping of all the generators and loads in the system. Moreover, it is assumed that all frequency-related control loops of the turbine-generators have the same structures and parameters. In this figure, Δ<sup>L</sup> is the changes in the loads, Δ<sup>m</sup> is the variations in the generator mechanical input, Δ<sup>c</sup> is the control power demanded by the secondary frequency control, Δ is the variation in the grid frequency, R is the droop coefficient, T<sup>1</sup> , T2 , T<sup>3</sup> and T<sup>4</sup> are the turbine and governor time constants. The inertia emulation using an inverter-based system is added to this block diagram with the emulated inertia constant of H, and the converter represented by a simple time constant of T<sup>c</sup> [110]. This time constant represents the delay and the bandwidth of the inverter. The parameters used for simulating this simplified model are given in Table 7.1.

Figure 7.2: Simplified block diagram of frequency control loops in a simplified power system with the addition of emulated inertia (model according to [267]).

Table 7.1: The parameters used for the simplified single-machine power system.


The response of the simplified power system model of Figure 7.2 to a small step increase in the load (ΔL) is shown in Figure 7.3 for different values of the emulated inertia constant (H). As seen, the increase in the emulated inertia reduces the maximum ROCOF and limits the frequency nadir, at the price of a slower recovery of the grid frequency. The increase in the frequency settling time, after recovering from the frequency nadir, is considered the main drawback of increasing the system inertia using inertia emulation techniques.

Figure 7.3: The effect of increasing the emulated inertia constant using converter-interfaced systems on (a) the frequency, and (b) the ROCOF during an under-frequency event.

Figure 7.4: The pole-zero map of the simplified power system with the increasing value of the emulated inertia constant (H).

This can also be observed from the pole-zero map of the same system for different values of inertia, shown in Figure 7.4. With the increase in inertia, the dominant system poles move towards the origin, which indicate a slower decaying response and a higher settling time. The pole-zero map has been drawn by neglecting the inverter time constant (T<sup>c</sup> ), as it is much smaller than the time constants of the turbine and governor. However, unlike the physical inertia of synchronous generators, emulated inertia is a control parameter which can be altered in real-time according to the system requirements. Therefore, to solve the problem of slow frequency recovery, adaptive inertia controllers have been proposed [265], [268]–[271], where the emulated inertia is reduced to near zero after the frequency nadir is reached for a faster recovery of the frequency . These adaptive inertia emulation techniques are later discussed in detail in section 7.2.2.

It can be seen in Figure 7.3, that the inertia has no influence on the steady-state frequency value. The steady-state frequency depends on the system damping (D) and the droop (R). A similar analysis on the effect of increasing the damping, rather than the inertia constant, shows that increasing the damping improves the frequency nadir and reduces the steady-state frequency error. The frequency nadir is determined by both, the system inertia and the damping, with the damping having a greater impact [265]. It can be concluded, that increasing the inertia and the damping can reduce the ROCOF and frequency deviation during imbalances of generation and demand. However, the inertia constant should be reduced when the frequency is recovering.

# **7.1.2 Flywheels as the Active Power Source for Inertia Emulation**

Initially proposed for wind turbines [272], the topic of providing inertia support using converter-interfaced systems has gained a great attention over the past decade. Significant contributions have been made by researchers working in this field, for instance in [273]– [281], where a variety of techniques and controllers has been proposed to emulate the inertia response of the synchronous generator using converter-interfaced systems during frequency deviations. Comprehensive reviews of many of these techniques can be found in [282] and [283].

Unlike the synchronous generators, where the kinetic energy stored in their large rotating rotors is used for providing the inertia response, converters do not inherently have a source of energy. Therefore, an active power source is required for providing the power demanded by the inertia emulation controller. In many of the proposed control designs for inertia emulation, for instance in [273]–[276], the active power source behind the inverter are simplified, assuming an ideal DC-link with an unlimited energy content. However, in practice, the dynamics and limitations of the active power source can also influence the inertia response of such systems. When the active power source is included, wind turbines [278], [279], PV systems [280], [281], DC-link capacitors [229], [284], and batteries [285], [286] are suggested to be the source of energy for inertia emulation. However, as summarized in Table 7.2, each of these solutions have their own limitations in providing the power required for emulating the inertia response. These limitations can change the response of the system from an ideal inertia response. In wind turbines, the speed recovery required after the first seconds of providing the inertia response, alters the response of the wind turbine from an ideal inertia response [272], [287]. A poorly conducted speed recovery can even cause a second frequency drop. Furthermore, inertia emulation may cause rotor stall or even trigger instabilities if the aerodynamics of wind turbines are not fully studied. Photovoltaic systems cannot participate in under-frequency events when being operated at their maximum power point. Therefore, additional storage systems or operating the system below the maximum power point is required, which is not ideal from an economical aspect. In an attempt to solve this, authors have also suggested to use the DC-link capacitors of the inverters as the active power source for inertia support [229], [284]. However, considering the low capacity of these capacitors and the voltage constrains for maintaining a linear pulse width modulation, the amount of energy that it can be provided for inertia emulation using these capacitors is very limited. Battery energy storage systems can be a good alternative for inertia emulation [285], [286]. However, concerns over the lifetime of the battery may force the operator to limit its power ramp rate, which can impair the inertia response. Most adaptive inertia emulation techniques, like the methods proposed in [264], [265], [268]–[271], require a sudden burst of power for several seconds, which is not the ideal use of a battery, concerning its lifetime and capacity. Therefore, it has been suggested to combine the battery with a supercapacitor in a hybrid structure, where the supercapacitor is mainly responsible for the inertia support [288]. However, this design leads to a higher degree of complexity in both control and hardware requirements. Moreover, the emulated inertia is limited by the size of the supercapacitor.

As discussed in Chapter 2, a FESS can rapidly provide high amounts of power with a high ramp rate and no concern over its lifetime or capacity. The FESS does not have any of the aforementioned limitations and can provide an ideal inertia response. In a highspeed FESS, the high-inertia rotor can store a significant amount of kinetic energy. However, the high-inertia rotor is electrically decoupled from the grid by a back-to-back converter. Nonetheless, by implementing an inertia emulation controller on the converter, the inertia can be transferred artificially to the grid side, providing an inertia response, similar to the response of the synchronous generators. The inertia emulation controller for a FESS is described in detail in section 7.2.


Table 7.2: Technical limitations of different active power sources for inertia emulation.

# **7.1.3 Inertia Emulation versus Virtual Synchronous Machines**

There are two common control structures for providing inertia support using converterinterfaced systems. In one category, the virtual synchronous machines [289]–[291], a low- or high-order model of a synchronous generator is used to generate the reference values for the voltage amplitude and phase angle of a grid-connected converter. In this design, the active power is controlled by changing the voltage phase angle, while the reactive power is controlled by varying the voltage amplitude, similar to synchronous generators. The main advantage of this approach is that, in some designs, it can eliminate the need for a frequency measurement system such as a phase-locked loop [291]. It uses the phase angle from the generator equation for synchronization of the converter. However, implementing virtual synchronous machines requires a drastic change and redesign of the converter controllers. Therefore, it is not logical or practical to implement such controllers on existing converters with already built-in controllers [279]. Also, such controllers do not inherently have current controllers. Therefore, an additional protection system to limit the currents may be required [292]. Moreover, for providing the inertia support using virtual synchronous machines, the active power controller including the inertia controller should be implemented on the grid-side inverter. However, in most commercial FESS, similar to many wind turbines, the active power controller is implemented on the machine-side converter. This design choice allows to directly control the machine torque, and facilitates parallel connection of several flywheels to reach higher ratings.

In the second category, which is referred to as inertia emulation techniques, only the swing equation of the synchronous generator is used to generate the active power reference for the converter. Therefore, the inertia emulation controller can be easily added as an outer control layer to the conventional current-controlled voltage source converters, including existing systems [292]. Due to this major advantage, this control design has been used by the European VSYNC project [293], as well as in [273], [279], [294], [295]. This method is also applied in this work, as it aims to implement the inertia emulation control as an external control loop to a commercial high-speed FESS, with preexisting converter controllers.

The drawback of this approach is the need for estimating the frequency and its derivative. In this work, the challenges associated with the use of a conventional phase-locked loop and the differential operator for estimating the frequency derivative are avoided by using a Dual Second-order Generalized Integrator Frequency-locked Loop (DSOGI-FLL). This is explained later in subsection 7.2.4.

### **7.1.4 State-of-the-Art in using a FESS for Inertia Emulation**

The use of a FESS for inertia emulation has been previously reported in [45] only, while the work in [46] presents a virtual synchronous generator using a FESS. However, in both cases, the parameters of the proposed controllers for inertia support are constant and adaptive controllers are not applied. The controller in [46], as a virtual synchronous machine controller, has no current control and authors report problems in forming pure sinusoidal waveforms. In both studies, realistic grid scenarios are not provided, and the performance of the proposed controllers has been demonstrated using numerical simulations only. To our knowledge, no experimental validation of using a FESS for inertia emulation has been reported. Furthermore, using PHIL testing for the validation of inertia emulation techniques has been reported only in [273], but the tests include only testing the inverter with a simple control design, and an ideal DC voltage source has been used as the active power source. Therefore, PHIL testing of a full system, including the active power source for inertia emulation, has not been conducted before.

This work attempts to improve the current state-of-the-art on using a FESS for inertia emulation by introducing a new adaptive inertia emulation controller, implementing it on a real FESS, and validating it using PHIL testing. The proposed controller combines the advantages of bang-bang controllers with self-adaptive ones for an improved response to frequency disturbances. Moreover, it considers the state of charge of the FESS in providing the inertia response and considers the practical limitations of using a FESS for frequency support. The proposed controller is described in the following section.

# **7.2 Adaptive Inertia Emulation Control for a Highspeed FESS**

In this section, the proposed adaptive inertia emulation controller for high-speed FESS is described in detail. Firstly, the controller structure for inertia emulation using FESS is explained, followed by the description of the suggested adaptive control method to alter the inertia constant and damping coefficient for an improved frequency regulation. Next, the choice of the parameters for the proposed controller is explained, considering the practical limitations of a FESS. Lastly, the method used to estimate the frequency and its derivative, as the inputs of the proposed inertia emulation controller, is described briefly.

### **7.2.1 Controller Design of Inertia Emulation using Highspeed FESS**

From the mechanical model of a PMSM in a high-speed FESS, described in Chapter 5, the electrical power of the PMSM (<sup>e</sup> ) is calculated as

$$P\_{\mathbf{e}} = \mathbf{J}\_{\mathbf{f}} \boldsymbol{\omega}\_{\mathbf{m}} \frac{d\boldsymbol{\omega}\_{\mathbf{m}}}{dt} + \mathbf{D}\_{\mathbf{f}} \boldsymbol{\omega}\_{\mathbf{m}}.\tag{7.1}$$

Equation (7.1) can be rewritten in per unit values using the definition of the inertia constant for a PMSM as [231]

$$\mathbf{H} = \frac{\frac{1}{\mathcal{Z}} \mathbf{J} \mathbf{\hat{o}}\_{\text{max}}^2}{\mathcal{S}\_{\text{n}}}.\tag{7.2}$$

In Eq. (7.2), max is the maximum speed of the FESS, and S<sup>n</sup> is its nominal apparent power. The inertia constant H given in seconds determines the time the flywheel can provide the nominal power using the stored kinetic energy only. By indicating the per unit value with an upper bar, Eq. (7.1) becomes

$$\overline{P}\_{\mathbf{e}} = 2\mathcal{H}\frac{d\overline{\omega}\_{\mathbf{m}}}{dt} + \mathcal{D}\overline{\omega}\_{\mathbf{m}}.\tag{7.3}$$

where D is the damping coefficient in per unit. Assuming the number of pole pairs of the PMSM to be one, the mechanical and electrical angular velocity become identical, therefore,

$$
\overline{P}\_{\mathbf{e}} = 2\mathcal{H}\frac{d\overline{\omega}\_{\mathbf{e}}}{dt} + \mathcal{D}\overline{\omega}\_{\mathbf{e}}.\tag{7.4}
$$

It is evident from Eq. (7.4) that the PMSM in a FESS is governed by the same swing equation of the conventional synchronous generators, with the difference of having no mechanical power input in a FESS. With any changes in the angular velocity of the supplied voltage to the PMSM, the PMSM provides an inertia response according to Eq. (7.4). However, in a FESS, the PMSM is decoupled from the grid by a back-to-back converter. To go around this problem, we replace the electrical angular velocity of the PMSM ̅<sup>e</sup> with the error in the grid's angular velocity in per unit Δ̅<sup>g</sup> . The resulting equation is then used to generate the active power reference (̅ ref) for the FESS. That is,

$$\mathcal{P}\_{\text{ref}} = 2\mathcal{H}\frac{d\Delta\overline{\alpha}\_{\text{g}}}{dt} + \mathcal{D}\Delta\overline{\alpha}\_{\text{g}}.\tag{7.5}$$

Equation (7.5) is the basis for the proposed inertia emulation controller. Using this approach, the inertia response of the PMSM is artificially transferred to the grid side, and the FESS reacts to the grid frequency deviations, similar to a synchronous generator.

The block diagram of the proposed inertia emulation controller is shown in Figure 7.5. As mentioned before, due to the structural design of commercial FESS, the active power control is often implemented on the MSC controller. The active power reference generated by Eq. (7.5) is first transformed to the electrical torque reference eref using the instantaneous speed of the flywheel m, which is then converted to the reference value for the q-axis current of the MSC using

$$\mu\_{\rm qref} = \frac{2\pi\_{\rm eref}}{3n\_{\rm p}\psi\_{\rm PM}} = \frac{2\bar{P}\_{\rm ref}\mathbf{P}\_{\rm n}}{3n\_{\rm p}\psi\_{\rm PM}\omega\_{\rm m}}.\tag{7.6}$$

The reference value for the d-axis current, the inner current controllers, and the GSC controller are identical to the ones presented in Chapter 5. The GSC controls the voltage of the DC-link according to the power injected or drawn from the MSC. Therefore, the power of the PMSM is transferred immediately to the grid. The GSC also regulates the voltage on the grid side, by controlling the reactive power according to the Q(U) characteristic of the German grid code, as shown in Figure 4.5(b).

Figure 7.5: The block diagram of the proposed inertia emulation controller, and the State of Charge (SOC) controller. This controller is to be added the MSC controller of the FESS in Chapter 5.

It can be observed from Figure 7.5, that the inertia constant H determines the amount of active power injection with respect to the ROCOF (Δ̅ <sup>g</sup> ), while the damping coefficient D determines the active power reference according to the error in frequency itself (Δ̅<sup>g</sup> ). As mentioned before, a major advantage of inertia emulation in comparison to the physical inertia is that the emulated inertia and damping are control parameters, which can be altered in real-time, for an improved response to frequency disturbances. As shown in subsection 7.1.1, the increase in the inertia constant can help to reduce the ROCOF and to improve the frequency nadir in low-inertia power systems. The effect of the damping coefficient resembles the role of the active power droop in a conventional frequency control system. The increase in the damping coefficient can improve the frequency nadir. In this work, a new adaptive control method is proposed for altering these two parameters in Eq. (7.5), with the aim of reducing the ROCOF and improving the frequency nadir in low-inertia power systems. Details of the proposed adaptive inertia and damping method are discussed in subsection 7.2.2.

It is also important to ensure that the FESS is always capable of providing and absorbing active power when a frequency disturbance occurs. Therefore, the FESS is also equipped with a State of Charge (SOC) controller, which attempts to maintain the state of charge of 75 %, after the frequency deviation has passed. The choice of 75% is based on the fact that flywheels are not able to provide the nominal power at low speeds, as previously discussed. It was demonstrated in subsection 6.3.1 that the nominal power can be delivered only above the SOC of charge of approximately 25 %. Therefore, 75 % is the center of the range, in which the FESS can provide its full power. The state of charge of 75 % corresponds to 86.6 % of the maximum speed of the FESS. Therefore, ref is set to 0.866max in Figure 7.5. The state of charge controller uses a conventional PI controller in combination with a ramp limiter, to limit the effect of charging or discharging of the FESS on the grid, and it is activated only when the frequency has recovered fully to the nominal value.

In the block diagram of Figure 7.5, a frequency deadband (±db) and a secondary frequency control block are included to provide a comprehensive frequency control framework. However, they are disabled in this work to show more clearly the effect of the inertia constant and the damping coefficient on the frequency response.

## **7.2.2 Adaptive Inertia and Damping for an Improved Response**

As mentioned before, an adaptive strategy can be used for the inertia constant and the damping coefficient to improve the response of the FESS to frequency deviations. In recent years, numerous techniques and solutions have been proposed in literature on how these parameters should be altered during a frequency disturbance. The proposed methods in literature can be categorized into two main groups: methods based on online optimization algorithms [265], [296]–[300], and interval-based controllers [268]–[271]. In this thesis, a new interval-based controller is proposed, and its performance is compared with several previously reported methods in this category. It is important to note that none of the previously reported adaptive methods are suggested specifically for a FESS. However, they can be applied to any converter-interfaced system. Therefore, they are used as reference points in this work for a comparison with the proposed method.

In interval-based controllers, the time period of the frequency deviation is divided into several fragments, according to the sign of Δ̅<sup>g</sup> Δ̅ g , as shown in Figure 7.6. When Δ̅<sup>g</sup> and Δ̅<sup>g</sup> have the same signs, i.e., Δ̅<sup>g</sup> Δ̅ g > 0, the period is called the accelerating phase. This includes the first interval after a disturbance occurs until reaching the frequency nadir. When they have opposite signs or Δ̅<sup>g</sup> Δ̅g ≤ 0, the interval is named the decelerating phase. In Figure 7.6, the accelerating phase is marked in red while the decelerating is colored in yellow. The naming of these periods resembles the behavior of a synchronous generators on its power-angle curve following a disturbance [268].

Figure 7.6: Segmentation of the frequency deviation period into the accelerating phase (red) and the decelerating phase (yellow), according to the sign of Δ̅<sup>g</sup> Δ̅ g .

As shown in subsection 7.1.1, having a large value for the inertia constant leads to a reduction of the maximum ROCOF observed during a frequency disturbance. However, maintaining the high inertia value leads to a slow recovery of the frequency, and an increased setting time for the frequency [268]. To solve the problem of the slow frequency recovery, a bang-bang control strategy for the inertia value is proposed in [268], where the inertia changes between two discrete values, depending on being in the accelerating or decelerating phases. Using this approach, when the frequency is moving away from the nominal value (accelerating phase), a high inertia is applied to resist its changes. As soon as the frequency begins to recover towards a new steady-state value (decelerating phase), the emulated inertia is set to zero to speed up the frequency recovery. However, the second degree of freedom, the damping coefficient is neglected in [268]. Therefore, in [269], a similar approach using the bang-bang control has been applied on both the damping and the inertia. The use of the damping in this design improves the frequency nadir and the frequency settling time. In both control designs using the bang-bang approach, the inertia constant and the damping coefficient can change only between two predetermined discrete values, which are independent from the severity of the frequency disturbance. In the adaptive inertia controller in [271], however, an adaptive component is added to the inertia, that scales with the amplitude of Δ̅<sup>g</sup> Δ̅ g . But this work, again, neglects the use of the damping coefficient. A self-adaptive controller is proposed in [270], which considers both the damping and the inertia constant. However, the damping is increased only after the frequency nadir (decelerating phase) to damp out the subsequent frequency oscillations. Therefore, it cannot influence the frequency nadir at all. The main issue with the proposed adaptive controllers in [271] and [270] is that the changes in the inertia and damping are not as fast as the bang-bang controllers, where the changes are almost immediate. In low-inertia power systems, where a high ROCOF is observed within the first instances of the frequency disturbance, it is important to reach a high inertia constant as fast as possible to effectively limit the maximum ROCOF. Enabling a high inertia value later than necessary can have little to no impact on the maximum ROCOF observed during the disturbance. A summary and a comparison of the interval-based adaptive controllers studied in this work is presented in Table 7.3.

literature. **Literature Bang-bang control Self-adaptive control Considering damping Considering limits of the active power source**

Table 7.3: A comparison of the proposed inertia emulation control with several interval-based methods in


\* The higher damping is applied only after the frequency nadir.

In this work, we combine the advantages of bang-bang control approaches, which guarantee a rapid change of the control parameters, with self-adaptive controllers, that ensure a more proportional response to the frequency disturbance. Moreover, we take advantage of the damping coefficient to help reducing the frequency nadir. The proposed adaptive control for the inertia constant and damping coefficient is described by Eq. (7.7) and (7.8), respectively.

$$H = \begin{cases} H\_1 + K\_H \left| \frac{d\Delta\overline{\alpha}\_\mathsf{g}}{dt} \right| & \text{, if } (\Delta\overline{\alpha}\_\mathsf{g} \frac{d\Delta\overline{\alpha}\_\mathsf{g}}{dt} > 0) \cap \left( \left| \frac{d\Delta\overline{\alpha}\_\mathsf{g}}{dt} \right| > \varepsilon\_\mathsf{H} \right) \\ & H\_2 & \text{, if } (\Delta\overline{\alpha}\_\mathsf{g} \frac{d\Delta\overline{\alpha}\_\mathsf{g}}{dt} \le 0) \cup \left( \left| \frac{d\Delta\overline{\alpha}\_\mathsf{g}}{dt} \right| \le \varepsilon\_\mathsf{H} \right) \end{cases} \tag{7.7}$$
 
$$D = \begin{cases} D\_1 + K\_D \left| \Delta\overline{\alpha}\_\mathsf{g} \right| & \text{, if } (\Delta\overline{\alpha}\_\mathsf{g} \frac{d\Delta\overline{\alpha}\_\mathsf{g}}{dt} > 0) \cap \left( \left| \frac{d\Delta\overline{\alpha}\_\mathsf{g}}{dt} \right| > \varepsilon\_\mathsf{D} \right) \\ & D\_2 + K\_D \left| \Delta\overline{\alpha}\_\mathsf{g} \right| & \text{, if } (\Delta\overline{\alpha}\_\mathsf{g} \frac{d\Delta\overline{\alpha}\_\mathsf{g}}{dt} \le 0) \cup \left( \left| \frac{d\Delta\overline{\alpha}\_\mathsf{g}}{dt} \right| \le \varepsilon\_\mathsf{D} \right) \end{cases} \tag{7.8}$$

As seen from Eq. (7.7), as soon the frequency disturbance occurs (accelerating phase), the inertia constant is instantaneously increased from a lower value (<sup>2</sup> ) to a higher value (<sup>1</sup> ), which is further increased by a proportional term ( | Δ̅g <sup>|</sup>), that scales linearly with the absolute value of the ROCOF. Therefore, using the proposed approach, the inertia constant increases rapidly using the bang-bang approach, and according to the severity of the frequency disturbance, an even higher inertia constant is applied. When the frequency begins to recover (decelerating phase), the inertia constant is reduced immediately to a value near zero (<sup>2</sup> ) to avoid the increase in frequency settling time. The condition to apply the higher inertia is being in the accelerating phase and not being in steady-state conditions, which is determined by | Δ̅ g <sup>|</sup> <sup>&</sup>gt; . The addition of the second condition, avoids unnecessary switching between the two states, which are observed in [268], where is set to zero. A low inertia value is also applied when the frequency has reached a steady-state condition, i.e., | Δ̅ g <sup>|</sup> <sup>≤</sup> <sup>ε</sup>H.

A similar controller design is suggested for the damping coefficient. However, the adaptive component of the damping coefficient (Δ̅<sup>g</sup> ) scales with the frequency error, rather than the ROCOF. Also, the adaptive term is maintained when the frequency is recovering, while the constant term is reduced from <sup>1</sup> to <sup>2</sup> . This helps improving the damping of the subsequent frequency oscillations after the nadir, while the reduction of the constant term helps to limit the energy drawn from the FESS.

The differences between the proposed controller and the reviewed adaptive methods in literature is summarized in Table 7.3. As seen in Table 7.3, the proposed controller also considers the limitations of the active power source used for inertia emulation, which is discussed in the next subsection.

### **7.2.3 Parameters for the Inertia Emulation Control using a FESS**

Many related works on adaptive inertia emulation techniques do not consider the limitations of the active power source behind the converter, and assume an ideal DC-link with unlimited energy resources, as in [273]–[276]. This section proposes a simple methodology to include the state of charge and the practical limitations of the FESS at low rotational speeds into the control design by adjusting the parameters of the proposed inertia emulation controller according to the state of charge of the FESS.

In the proposed control design for the inertia constant, shown in Eq. (7.7), <sup>1</sup> and are the control parameters that determine the emulated inertia, which are used to calculate the active power injection according to the ROCOF. When the FESS is fully charged at its maximum speed of max, these two parameters should be at their maximum level to provide the maximum inertia. However, it is proposed in this work that these parameters are scaled down with the reduction of the state of charge of the FESS. Using this approach, the FESS will provide inertia support proportional to its state of charge.

Therefore, the inertia constant <sup>1</sup> is equal to its maximum value for a fully charged FESS, H1max, and for lower values of state of charge,

$$H\_1 = \frac{\frac{1}{2}\mu\omega\_\mathrm{m}^2}{\frac{1}{2}\mu\omega\_\mathrm{max}^2}H\_{1\max} = (\frac{\omega\_\mathrm{m}}{\omega\_\mathrm{max}})^2 H\_{1\max} = \mathrm{SOC} \cdot H\_{1\max}.\tag{7.9}$$

Figure 7.7(a) shows the inertia constant <sup>1</sup> as a function of the speed of the FESS according to Eq. (7.9), and for a maximum inertia constant of 5.9 s. A similar function to Eq. (7.9) is used for with the maximum value of Hmax, i.e., = SOC ∙ Hmax. The parameter H<sup>2</sup> is constant and close to zero to speed up the frequency recovery.

Furthermore, we integrate the power limitations of the FESS at low rotational speed into the controller. This limitation is integrated into the damping coefficient <sup>2</sup> , as it determines the FESS contribution in steady-state conditions. The damping coefficient <sup>2</sup> is reduced from its maximum value 2max when the FESS reaches the low speed threshold, according to the speed of the FESS, as shown in Figure 7.7(b). A similar function is used for <sup>1</sup> and , with the maximum values of 1max and Dmax, respectively, as shown in Figure 7.7.

Using the proposed approach, the FESS contribution to frequency deviations is reduced automatically, when the FESS reaches lower speeds and reaches zero, when the FESS is emptied. This is a mandatory requirement, when using a real FESS, where the physical system constrains must be considered.

As the damping coefficient <sup>2</sup> determines the steady-state contribution of the FESS, it has the same effect as the droop coefficient. Therefore, we chose the maximum value of 2 to be 40 p.u., which corresponds to a droop value of 2.5 %. The value for <sup>1</sup> is chosen to be slightly higher to help to improve the frequency nadir, as discussed earlier. The choice of Hmax and Dmax determines the effect of the instantaneous ROCOF and frequency on the inertia constant and the damping coefficient, respectively. These are design parameters, which can be selected according the maximum ROCOF and frequency deviation expected in a system. The value of ε<sup>H</sup> and εD, which determine the transition condition in the proposed controller, should be higher than the noise level expected on the ROCOF and frequency estimation, in order to avoid unwanted triggering between the two control states. In this respect, using a differential operator to calculate the ROCOF can lead to noise amplification and inaccurate triggering of the controller. Therefore, in the next section, a method is described to calculate the ROCOF without using this operator. The parameters used for the proposed inertia emulation control in this work are given in Table 7.4.

Figure 7.7: The variations of the inertia emulation control parameters with respect to the speed of the FESS. (a) The inertia constants, <sup>1</sup> and <sup>2</sup> , (b) the damping coefficients <sup>1</sup> and <sup>2</sup> , and (c) the adaptive coefficients <sup>H</sup> and D.


Table 7.4: The parameters used for the proposed inertia emulation controller.

#### **7.2.4 Estimation of the Frequency and its Derivative**

Inertia emulation techniques require estimated values of the frequency and the frequency derivative (ROCOF). The Synchronous Reference Frame Phase-locked Loop (SRF-PLL) has been used often in inertia emulation controllers in literature [290], [294], [295], [300], as the conventional method for estimating the frequency. However, it has been reported that the use of the SRF-PLL in inertia emulation applications can lead to instabilities and noise amplification associated with the differential operation, used for calculating the ROCOF [301]–[303]. Furthermore, the SRF-PLL performs poorly under unbalanced and distorted voltage conditions, often observed in weak grids, such as microgrids [304]. Although the performance of the SRF-PLL in such conditions can be improved by additional filters or moving window functions, these methods introduce more delays into the frequency estimation. For the application of inertia emulation in low-inertia power systems, a fast detection of frequency changes is required, and such delays can impair the benefits of inertia emulation.

Therefore, in this work, it is suggested to use the Dual Second-Order Generalized Integrator Frequency-Locked Loop (DSOGI-FLL), introduced in [304]. Using the DSOGI-FLL has the significant advantage that it can estimate the frequency derivative without using the differential operator, avoiding the problems of noise amplification and unrealistic outputs. The DSOGI-FLL only uses integrators and multiplication of filtered signals to estimate the ROCOF. The superior performance of the DSOGI-FLL over the SRF-PLL for the estimation of the ROCOF has been demonstrated by simulations and experiments in [302]. The same work shows that the common stability issues of the SRF-PLL are not observed in the DSOGI-FLL. In addition, the DSOGI-FLL is based on adaptive filters. Therefore, it performs well under distorted, faulty, and unbalanced voltage conditions [304]. The block diagram of the DSOGI-FLL is depicted in Figure 7.8.

Figure 7.8: The block diagram of the Dual Second-Order Generalized Integrator Frequency-Locked Loop (DSOGI-FLL) (adopted from [304]).

As seen in Figure 7.8, the DSOGI-FLL consists of two Second-order Generalized Integrator Quadrature Signal Generators (SOGI-QSG) and a Frequency-locked Loop (FLL), with a gain normalization block. The two SOGI-QSG receive the measured voltage in the stationary reference frame ( and ), and generate filtered in-phase (′ and ′) and in-quadrature (′ and ′) signals of the input voltages with the same amplitude. The transfer function of these conversions for the in-phase ( ()) and for the in-quadrature ( ()) components are shown in Eq. (7.10) and (7.11), respectively, which are identical for both and components.

$$G\_{a\emptyset}^{d}\{\mathbf{s}\} = \frac{\stackrel{\cdot}{v}\_{a\emptyset}}{v\_{a\emptyset}} = \frac{\text{k}\omega\_{\emptyset}s}{\text{s}^2 + \text{k}\omega\_{\emptyset}s + \omega\_{\emptyset}^2}.\tag{7.10}$$

$$G\_{a\theta}^{q}(\mathbf{s}) = \frac{q\nu\_{a\theta}^{\prime}}{\nu\_{a\theta}} = \frac{\mathbf{k}\omega\_{\mathbf{g}}^{2}}{\mathbf{s}^{2} + \mathbf{k}\omega\_{\mathbf{g}}\mathbf{s} + \omega\_{\mathbf{g}}^{2}}.\tag{7.11}$$

In these equations, <sup>g</sup> is the estimated frequency by the FLL and k is the SOGI-QSG gain, which is a design parameter. By looking at Eq. (7.10) and (7.11), it can be seen that the SOGI-QSG control loops are basically a band-pass filter for the in-phase components and a low-pass filter for the in-quadrature signals. To have the optimum damping ratio () of 0.707, the SOGI-QSG gain k is set to √2, as recommended in [304]. This choice leads to an optimum balance between a fast dynamic response and the filtering capability of the SOGI-QSG. The FLL uses the error in calculating the in-phase components and the filtered in-quadrature components to calculate the frequency derivative as

$$\frac{d\overline{u}\_{\overline{\mathbf{g}}}}{dt} = -\frac{\Gamma}{2}(\{\boldsymbol{\nu}^{\prime}\boldsymbol{\boldsymbol{\nu}} - \boldsymbol{\nu}\_{a}\}q\boldsymbol{\nu}^{\prime}{}\_{a} + \{\boldsymbol{\nu}^{\prime}\_{\beta} - \boldsymbol{\nu}\_{\beta}\}q\boldsymbol{\nu}^{\prime}{}\_{\beta}).\tag{7.12}$$

In Eq. (7.12), Γ is the FLL gain, which is assumed to be equal to 100, as recommended in [251]. This choice of Γ leads to a relatively fast frequency setting time of the frequency measurement. It is clear from Eq. (7.12) that no differential operation is involved in estimating the ROCOF, when using the DSOGI-FLL. The gain normalization block ensures that FLL receive signals with unity magnitude, avoiding errora in frequency estimation [251]. The frequency itself is estimated by integrating over the calculated frequency derivative.

# **7.3 Simulation Results of the Adaptive Inertia Emulation using a high-speed FESS**

The performance of the proposed adaptive inertia emulation controller for the high-speed FESS is initially evaluated using offline simulations in MATLAB/Simulink environment. Two difference scenarios in a low voltage microgrid are presented to show the effectiveness of the proposed method in improving the frequency dynamics during power imbalances. For a comparison with the current state-of-the-art, three of the interval-based adaptive inertia controllers, previously reported in the literature, are also simulated in each scenario.

The simulation scenarios are carried out in a low voltage microgrid, based on the CIGRE European low voltage network benchmark [123], and its variation in [305]. Among the available benchmark grids with openly available data, this network benchmark has the closest characteristics to a typical German low voltage distribution grid. The high-speed FESS together with two PV systems and a diesel-based Combined Heat and Power (CHP) system are added to the CIGRE benchmark to form a microgrid and allow an autonomous operation of the system. The single-line diagram of the microgrid is shown in Figure 7.9. The FESS is installed near the transformer to provide frequency and voltage support during islanding, while the PV systems and the CHP are located near the loads. The microgrid is capable of being disconnected from the main supply and operate in the autonomous mode, when necessary. By a power imbalance in the autonomous operation of the microgrid, the frequency can deviate greatly, as the synchronous generation of the CHP is the only source of traditional inertia in the system. During an imbalance of generation and demand, the CHP slowly covers the demand, while the FESS has the task of providing a faster response and additional inertia support during the first seconds of the frequency disturbance. The PV systems are modelled according to the generic two-stage model presented in [261] and are being operated at their maximum power point. Therefore, they are not able to take part in regulating the frequency during under-frequency events. The parameters of the diesel-based CHP and the PV systems are given in the Appendix B. The FESS parameters are the same as the ones use in Chapter 5 and 6, and the inertia emulation parameters are given in Table 7.4.

Figure 7.9: The single-line diagram of the test microgrid based on the CIGRE European low voltage network benchmark.

Two different scenarios are simulated, which are the microgrid islanding and the intermittent generation of the PV systems in the autonomous mode of operation. In each scenario, six difference cases are simulated:


A summary of the characteristics of the controllers in all cases is presented in Table 7.5. For a fair comparison, the same maximum inertia constant and damping of Table 7.4 is used for the bang-bang control in case 3. For the controllers in case 4 and 5, where a single value of inertia constant is used, the average of inertia and damping values of Table 7.4 is used. The simulation results in the two scenarios and for all the six cases are provided next.


Table 7.5: A summary of the controller characteristics in all cases.

#### **7.3.1 Scenario A: Microgrid Islanding**

In this scenario, the microgrid is initially drawing 80 kW of power from the medium voltage grid, while the PV systems and the CHP system are only partially supporting the loads. The FESS is charged to reach the state of charge of 75 % using the state of charge controller, prior to the disturbance. The microgrid is decoupled instantaneously from the main supply by opening the circuit breaker S1, and it undergoes into the islanded mode of operation. The power deficiency in the microgrid leads to a sudden drop in the microgrid frequency. The FESS attempts to limit the ROCOF and frequency nadir, while the CHP system slowly increases its output to compensate for the power deficiency.

The simulation results for all the six cases in this scenario are presented in Figures 7.10 and 7.11. It can be seen from Figures 7.10(a) and (b), that the use of the FESS with the proposed adaptive inertia emulation controller, presented as case 6, leads to the highest reduction in the maximum ROCOF and the most improved frequency nadir. As shown in Figure 7.10(b), the FESS with the proposed controller reduces the maximum ROCOF by 28 %, from 1.72 Hz/s in case 1, to 1.24 Hz/s in case 6. Furthermore, it performs better than the previously proposed adaptive inertia controllers, presented as case 3 to 5, which reduce the maximum ROCOF to 1.34 Hz/s at best. The superior performance of the proposed design in reducing the ROCOF is due to quickly reaching a higher inertia constant, as shown in Figure 7.10(c), using a combination of the bang-bang control approaches and the self-adaptive ones. In case 3, where only the bang-bang control is used, the inertia constant also changes very rapidly. However, it only reaches a predetermined constant value (1max), independent from the ROCOF, and remains constant until the frequency nadir is reached, which limits its effectiveness. The adaptive inertia controller implemented in case 4 changes the inertia constant slowly, and therefore, it performs worse than the proposed controller. In the self-adaptive controller of case 5, the inertia constant is only increased by the instantaneous value of the ROCOF. Therefore, the increase is not as high and as fast as in the proposed design. As a result of combing the bang-bang control and the self-adaptive approach in the proposed design, the inertia constant immediately reaches the high value of 4.42 s (<sup>1</sup> at SOC of 75%), which is increased even further by the instantaneous value of the ROCOF, reaching an inertia constant of 7.8 s quickly after islanding. Among the previously reported control designs, case 3 and 5 show a similar behavior in reducing the maximum ROCOF. The controller in case 5 reduces the ROCOF slightly better than the one in case 3, due the fact that the initial ROCOF in this scenario is high enough to reach a higher inertia constant than the pre-given value in case 3. Although a higher inertia constant is reached in case 4, in comparison to case 3, the increase is not fast enough to effectively limit the maximum ROCOF. Thus, the worst performance is observed in case 4 among the cases 3 to 5. The maximum ROCOF and frequency nadir in all cases are summarized in Table 7.6.

The proposed controller also results in the largest improvements in the frequency nadir, by taking advantage of the second degree of the freedom, the damping coefficient. As seen in Figure 7.10(a), the use of the FESS with a conventional droop control, presented as case 2, improves the frequency nadir from 48.84 Hz, to 49.2 Hz. However, with a proper use of the damping coefficient in the proposed design, the frequency nadir can be improved even further, to reach 49.35 Hz. This is a 44 % reduction in the maximum frequency deviation in comparison to case 1, and 18.7 % in comparison to case 2. This is achieved by increasing the damping coefficient immediately after the islanding, and further increasing it according to the instantaneous changes in the frequency, as shown in Figure 7.10(d). Between the studied cases, the controller in case 4 does not use the damping coefficient, leading do the worst frequency nadir. Following case 6, the most improved frequency nadir can be observed in case 3, as unlike the controller in case 5, a higher damping coefficient is applied before reaching the frequency nadir. It should also be noted, that the improved frequency nadir is also a side effect of the reduction in the ROCOF in the proposed design using the inertia constant. Using the FESS with the proposed control design also shows a better damping of frequency oscillations after the frequency nadir, as seen in Figure 7.10(a). This is due to maintaining the adaptive component of the damping coefficient, after reaching the frequency nadir.

Figure 7.10: Simulation results for scenario A: microgrid islanding. (a) Frequency, (b) ROCOF, (c) Inertia constant (H), and (d) Damping coefficient (D).

Figure 7.11 shows the active power of the FESS and its state of charge variations during the first five seconds of the islanding scenario. As seen in Figure 7.11(a), the use of the proposed controller leads to a more prompt and faster response of the FESS, which leads to the improved grid frequency dynamics during the islanding. The controller design in case 3 results in a slightly higher peak in power than the proposed design; however, a

better performance in terms of the frequency is not observed, as the increase in power is slower than in the proposed control design. In case 3 and 6, where the bang-bang control and a higher damping coefficient from the beginning of the frequency incident are applied, a higher peak in the FESS power is observed. However, the power is reduced significantly using the bang-bang approach, after reaching the frequency nadir. Therefore, the total energy drawn from the FESS in all cases do not differ significantly, as seen in Figure 7.11(b). The proposed control design increases the energy drawn from the FESS by only 0.15 %, in comparison to case 2, and by 0.07 %, in comparison to case 3 and 5. It can be concluded, that the proposed adaptive inertia emulation controller leads to an improved frequency dynamic of the microgrid, in terms of reducing the maximum ROCOF and improving the frequency nadir, while not using significantly more energy from the FESS. This is an advantage in a long-term autonomous operation of the microgrid with continuous variations in load and generation.

Figure 7.11: Simulation results for scenario A: microgrid islanding. (a) FESS Power, and (b) FESS State of Charge (SOC).

Table 7.6: The maximum ROCOF and frequency nadir for all the simulated scenarios and cases.


#### **7.3.2 Scenario B: Drop in the PV Generation**

In scenario B, the microgrid is disconnected from the main supply, after an islanding scenario, similar to the one in scenario A. The CHP system and the two PV systems are supporting the loads, while the FESS is injecting a small amount of power according to the value of <sup>2</sup> and the steady-state frequency error. A sudden sharp drop in the solar irradiance density is simulated, reducing it from 1000 W/m<sup>2</sup> to 150 W/m<sup>2</sup> , which could be caused by passing clouds. This results in a drastic drop of the PV generation in the microgrid, from 70 kW to 10.5 kW. The power deficit leads to a drop in the frequency of the microgrid. The FESS increases rapidly its output to limit the ROCOF and frequency nadir, while the CHP system slowly covers the lack of generation by the PV systems.

The simulation results for the six cases in this scenario are shown in Figures 7.12 and 7.13. Similar to scenario A, the use of the FESS with the proposed controller leads to the highest reduction in the maximum ROCOF and frequency nadir during this frequency disturbance. It reduces the maximum ROCOF from 1.15 Hz/s in case 1, to 0.82 Hz/s in case 6. It can be seen from Figure 7.12(c) that this is achieved by quickly reaching a higher inertia constant in the proposed controller in comparison to the other cases. In contrast to the scenario A, the controller in case 3 performs better than case 5 in this scenario. This is due to the fact, that the inertia constant in case 5 increases only with the ROCOF, and since the ROCOF in this scenario is lower than the one in scenario A, case 3 with a predetermined increase in the inertia constant shows a better performance. Despite eventually reaching a higher inertia constant in case 5, the increase is not fast enough to limit the maximum ROCOF.

Similarly, the highest improvement in the frequency nadir is seen using the proposed control design, where the frequency nadir reaches only 49.34 Hz, in comparison to 49.04 Hz in case 1, and 49.22 Hz in case 2. Following case 6, case 3 shows the best performance in terms of the frequency nadir, by reaching 49.29 Hz. This is due to the fact that it also uses a higher damping coefficient before reaching the frequency nadir.

Figures 7.13(a) and (b) show the active power and the state of charge of the FESS, respectively. As seen in Figure 7.13(a), the FESS is injecting about 12 kW of power, prior to the disturbance. Similar to the previous scenario, with the proposed controller, the FESS increases its power output faster than in the other cases. It also reaches a higher peak in power in this scenario, but due to the reduction of the power after reaching the frequency nadir, there is only a slight increase in the energy required from the FESS, when using the proposed controller. This can be seen by comparing the state of charge before the disturbance with the one after the five seconds into the disturbance. It is important to note that the state of charge values in different cases are different prior to the disturbance, as the drop in the PV generation is simulated following an islanding scenario using different controllers.

Figure 7.12: Simulation results for scenario B: drop in PV generation. (a) Frequency, (b) ROCOF, (c) Inertia constant (H), and (d) Damping coefficient (D).

Figure 7.13: Simulation results for scenario B: drop in PV generation. (a) FESS Power, and (b) FESS State of Charge (SOC).

# **7.4 PHIL Testing of a High-speed FESS for Inertia Emulation**

To experimentally validate the proposed inertia emulation control for the FESS, the controller has been implemented on the 60 kW high-speed FESS, using the concept of rapid control prototyping. The performance of the FESS using the proposed adaptive inertia emulation controller has been evaluated using PHIL testing. The two simulated scenarios, i.e., the microgrid islanding and the drop in PV generation, have been repeated in real-time for conducting the PHIL testing of the FESS. In each scenario, all six cases have been investigated, leading to 12 experiments in total.

#### **7.4.1 Real-time Simulation of the Microgrid**

For conducting the PHIL testing, a real-time simulation of the grid under study is mandatory. The CIGRE-based low voltage microgrid shown in Figure 7.9 is simulated in real-time on the Opal-RT's OP5700 real-time simulator with simulation time steps of only 24 µs. This is achieved by using the State-space Nodal (SSN) solver [109] to decouple the model at bus R4, creating three groups of equations that can be solved almost independently (see subsection 3.1 for more information), and neglecting the line capacitances due to their short lengths [117]. The relatively low simulation time step leads to a low total loop delay, which is a main contributor to the PHIL stability, as discussed in subsection 3.2.2.

## **7.4.2 Implementing the Inertia Emulation Controller on the real FESS**

The same PHIL setup and commercial high-speed FESS, which were described in Chapter 6, are also used for the experiments in this chapter.

To implement the proposed inertia emulation controller on the commercial FESS, the controller is also simulated in real-time on the OP5700 real-time simulator, as an external controller. It measures the frequency and frequency derivative from the grid voltage at bus R11, and generates the active power reference for the high-speed FESS. This concept, where the controller is simulated in real-time, while it is being tested on a real hardware is known as rapid control prototyping. The active power reference generated from the simulated controller is sent in real-time to the FESS internal controller through a 4-20 mA analog channel. In this work, the analog connection has been preferred to the standard Modbus communication, as the Modbus communication in known to be slow and nondeterministic in signal transmission, which is not suitable for this application. A schematic of the described signal connections and the PHIL testing setup is shown in Figure 7.14.

Figure 7.14: The signal connections and PHIL setup for implementing the inertia emulation controller on the FESS using the concept of rapid control prototyping.

In Figure 7.14, R11 is the three-phase grid voltage at bus R11, where the FESS is virtually connected. This voltage is sent via a dedicated SFP connection to the Egston 200 kVA power amplifier, which is feeding the FESS. It is also sent within the simulator to the DSOGI-FLL and the inertia emulation controller. The DOSGI-FLL receives this voltage and estimates the frequency and frequency derivative, as described in subsection 7.2.4. The estimated frequency and ROCOF are sent to the simulated inertia emulation controller, where the active power reference (ref) is calculated. The inertia emulation controller also receives the state of charge of the FESS form the FESS controller via the Ethernet connection and the Modbus protocol. The active power reference is sent from the real-time simulator to the FESS controller via an analog connection. The FESS controller sends the active power reference to the MSC controller, where it generates the torque reference from the active power reference. The FESS injects or absorbs power, accordingly to the received set point. The FESS currents (FESS) are measured at the FESS terminals, sent to the real-time simulator via another SFP connection, and inserted into the simulated grid at bus R11. Similar to the PHIL tests in Chapter 6, the voltage-type ideal transformer method [256] is used for interfacing the FESS and the real-time simulation, and a first-order low-pass filter with a cut-off frequency of 100 kHz is applied on the current measurements fed back to the real-time simulation.

To illustrate the functionally of the setup, steps in the active power reference are sent from the real-time simulator to the FESS and the response of the FESS to the received set point is measured. The results are shown in Figure 7.15. First, a 60 kW step in the active power reference is sent to the FESS, while initially the FESS is injecting no power. It can be observed from Figure 7.15(a) that the FESS reaches the active power set point of 60 kW in just over 25 ms, which indicates a very quick response from FESS. Next, while the FESS is being charged at its nominal power, the active power set point is changed stepwise from -60 kW to +60 kW, resulting in a 120 kW step change in the active power. As shown in Figure 7.15(b), it takes less than 32 ms for the FESS to go from charging at full power to discharging at the same rate.

Figure 7.15: Response of the FESS to steps in the active power reference, sent from the OP5700 real-time simulator. (a) A 60 kW step, from steady-state operation and no power injection. (b) A 120 kW step, from charging at the nominal power to discharging at the same rate.

The same scenarios and cases, described in the simulations of section 7.3, are repeated using the real-time simulation of the microgrid and the real FESS. The obtained PHIL testing results are presented in the next subsections.

### **7.4.3 Scenario A: Microgrid Islanding**

In this scenario, the simulated microgrid, running in real-time on the OP5700 real-time target, is decoupled from the main supply. The medium voltage grid is providing 80 kW of power to the microgrid, prior to the islanding. For all the six cases, before the islanding, the FESS is charged using the state of charge controller to reach the state of charge of 75 %, which corresponds to a rotational speed of approximately 39,000 rpm.

The results of the PHIL testing of the FESS for this scenario are shown in Figures 7.16 and 7.17. The frequency and the ROCOF of the microgrid during the islanding are depicted in Figure 7.16(a) and (b), respectively. It can be observed, that similar to the simulation results, the best performance in terms of reducing the maximum ROCOF and improving the frequency nadir is achieved using the proposed adaptive inertia emulation controller. The use of the high-speed FESS with the proposed inertia emulation controller, reduces the maximum ROCOF during the islanding scenario by 25 %, from 1.72 Hz/s in case 1, to 1.3 Hz/s in case 6. The superior performance of the proposed control design is achieved by rapidly reaching a higher inertia constant, as shown in Figure 7.16(c), by a combination of the bang-bang control methods and self-adaptive ones. The inclusion of the adaptive component in the inertia constant of the proposed controller leads to a better performance in comparison to the simple bang-bang control of case 3, where the inertia values in each interval are constant. Similarly, exploiting the bang-bang control approach in the proposed design results in a faster change of parameters in comparison to the self-adaptive controller of case 4.

Among the previously reported adaptive inertia controllers, case 3 and 5 show similar results, reducing the maximum ROCOF to only 1.4 Hz/s. Despite reaching a higher inertia constant in case 4, in comparison to case 3, a further improvement of the ROCOF is not observed, due to the gradual change of inertia constant in self-adaptive controllers in comparison to the bang-bang control methods.

Furthermore, as shown Figure 7.16(a), the FESS with the proposed controller improves the frequency nadir from 48.84 Hz in case 1, to 49.27 Hz in case 6. This is achieved by exploiting the damping coefficient immediately after the disturbance occurs, as opposed to case 5, and in combination with quickly reaching a higher inertia constant. As discussed in section 7.1.1, the frequency nadir is determined by both the damping and the inertia, with the damping having a greater influence. After the proposed controller, the controller implemented as case 3 shows the best performance in terms of the frequency nadir, as it is the only method which employs a higher damping coefficient, before reaching the frequency nadir. The highest ROCOF and the worst frequency nadir are observed in case 4, since it does not use the damping coefficient, and the inertia constant is increased slowly.

Moreover, it can be seen in Figure 7.16(a), that using the proposed approach, a better damping of the subsequent frequency oscillations is observed, having almost no overshoot in the frequency, while recovering.

Figure 7.17(a) shows the measured active power of the FESS. To calculate the active power, the FESS voltage and current measurements are collected at its terminal using voltage and current transformers. The instantaneous active power is calculated from the voltage and current measurements in the real-time simulator.

Figure 7.16: PHIL testing results for scenario A: microgrid islanding. (a) Frequency, (b) ROCOF, (c) Inertia constant (H), and (d) Damping coefficient (D).

As seen in Figure 7.17(a), similar to the simulation results, the proposed controller results in a more rapid and prompt response from the FESS to the frequency disturbance, which

leads to the improvements observed in the frequency dynamics. The fast response of the FESS is due to the quick change of the inertia constant and damping coefficients using the bang-bang control approach. The active power in case 3 and 6, where the bang-bang approach on both the inertia and the damping are used, reaches a higher peak value. However, as soon as the frequency nadir is reached, the active power is reduced by more than 10 kW and reaches values lower than the power of the FESS in other cases, including the droop controller of case 2. Therefore, the energy drawn from the FESS is not increased significantly. The reduction in the state of charge of the FESS during the five seconds of the islanding scenario is depicted in Figure 7.17(b). In comparison to the conventional droop control, presented as case 2, the proposed inertia controller draws only 0.12 % more energy from the FESS, and only 0.05 %, in comparison to case 3. Therefore, it can be concluded that the improved frequency control using the proposed controller is achieved by only a negligible increase in the energy required from the FESS.

In addition, the PHIL testing results confirm the simulation results in this scenario, regarding the performance of the proposed controller with minor differences in the results. Slightly higher ROCOF values and lower frequency nadirs are observed in the PHIL testing in comparison to the pure simulation studies (see subsection 7.4.5).

Figure 7.17: PHIL testing results for scenario A: microgrid islanding. (a) Measured FESS power, (b) FESS State of Charge (SOC).

#### **7.4.4 Scenario B: Drop in the PV Generation**

In this scenario, a sudden drop in the solar irradiance density is simulated, from 1000 W/m<sup>2</sup> to 150 W/m<sup>2</sup> , while the microgrid is disconnected from the mains. The drop in the solar irradiance density, as the main input of the two PV systems, leads to a drastic decrease in the generated power of these systems. With the drop in PV generation, the frequency falls, while the CHP system slowly increases its power to compensate for the power deficit. The FESS attempts to respond as fast as possible to reduce the maximum ROFOC during the disturbance and limit the frequency nadir.

The results of the PHIL testing of the FESS for the six cases in this scenario are shown in Figues 7.18 and 7.19. Similar results in comparison to the simulations are obtained. It can be seen from Figure 7.18(b) that the use of the FESS with an inertia emulation controller, independent of the adaptive inertia and damping method, is able to reduce the ROCOF to below 1 Hz/s. However, the highest reduction in the maximum ROCOF is achieved using the proposed adaptive controller, from 1.16 Hz/s in case 1, to 0.83 Hz/s in case 6. Between the previously reported methods, the controller in case 3 performs better than the others, reducing the ROCOF to 0.9 Hz/s. Unlike scenario A, the controller in case 5 performs worse than the controller in case 3, due to the fact that the severity of the frequency disturbance is smaller in this scenario, and the control parameters in case 5 scale only with the amplitude of the ROCOF and the frequency error.

The frequency variations of the microgrid during this scenario are shown in Figure 7.18(a). It can be observed that the proposed controller improves the frequency nadir from 49.04 Hz in case 1, to 49.31 Hz in case 6, achieving the best results in terms of the frequency nadir. Moreover, a better damping of the frequency oscillations is observed, when the frequency is recovering to its pre-disturbance value.

The measured active power of FESS and its state of charge are shown in Figure 7.19. In comparison to scenario A, the FESS is injecting 12 kW of power, prior to the disturbance, according to damping coefficient <sup>2</sup> . The initial state of charge of the FESS is different in each case, as this scenario has been simulated after an islanding scenario using different controllers in each case. Again, using the proposed adaptive inertia emulation controller, the FESS increases its power more quickly, as seen in Figure 7.19(a), which leads to the improved frequency dynamics of the microgrid. However, the increase in the energy required from the FESS by using the proposed controller is not significant.

It can be concluded that the PHIL testing results of the 60 kW high-speed FESS with the suggested adaptive inertia emulation control confirm the offline simulation results. The proposed controller for the FESS can effectively limit the ROCOF and frequency nadir during an imbalance between generation and demand. Moreover, it outperforms several previously reported methods in literature, while not requiring noticeably more energy from the FESS. The use of the suggested control design requires a higher ramp rate of the active power drawn from the FESS when a disturbance occurs, which is not a problem for a FESS as a storage technology, as it has no effect on its lifetime or capacity.

Figure 7.18: PHIL testing results for scenario B: drop in PV generation. (a) Frequency, (b) ROCOF, (c) Inertia constant (H), and (d) Damping coefficient (D).

Figure 7.19: PHIL testing results for scenario B: drop in PV generation. (a) Measured FESS power, (b) FESS State of Charge (SOC).

## **7.4.5 Comparison of the PHIL Testing and the Simulation Results**

By comparing the simulation results presented in section 7.3 and the results of the PHIL testing of the real FESS in section 7.4, a good match between the two can be observed. However, slight differences can be seen, which could have been caused by the delays and limitations in the real hardware, which are not considered in the simulated model. Some examples include the nonlinear torque ramp limiter of the permanent magnet synchronous machine and the delays in signal transmission, which exist only in the real hardware, and are difficult to implement in the simulations. These delays and limitations result in slightly higher ROCOF values, lower frequency nadir, and higher power drawn from the FESS. This is a clear advantage of the PHIL testing as a validation tool, where the behavior of a real hardware including all its limitations can be directly reflected on the grid variables, such as the frequency.

In order to see the exact differences between the response of the real hardware and the simulated model of the FESS, the same active power reference, sent to the real FESS during the PHIL testing, is sent simultaneously to the real-time simulation model of the FESS, developed in Chapter 5. The model is evaluated by comparing the measured active power of the real FESS and the simulated one. As an example, the results for the scenario B using the proposed controller (case 6) are shown in Figure 7.20. As seen, there is a good match between the response of the real-time model and the real FESS with a slight difference at the peak power. The difference between the peak power of the model and the real FESS is only 0.93 kW, which corresponds to only 1.8 % of error.

Figure 7.20: A comparison of the real FESS behavior and its simulated model, given the same active power reference from the PHIL testing.

# **7.5 Summary**

Low-inertia power systems can experience a high ROCOF and large frequency deviations during imbalances in generation and demand. In this chapter, a new adaptive inertia emulation controller for a high-speed FESS is introduced, aiming to simultaneously reduce the ROCOF and improve the frequency nadir in low-inertia systems during disturbances. The proposed adaptive method combines the advantages of a bang-bang controller and a self-adaptive one, in order to have a fast and yet proportional response to the frequency disturbances. The proposed control design also has the advantage that it can be easily added as an external control to a commercial FESS with no change required in the existing control systems. Furthermore, the proposed design integrates the state of charge and limitations of the FESS at low rotational speeds into the control architecture.

The performance of the proposed controller has been initially evaluated using numerical simulations in a low voltage microgrid. The controller is also implemented on a real 60 kW high-speed FESS using the concept of rapid control prototyping. The performance of the FESS with the proposed adaptive inertia emulation controller, along with several previously suggested controllers, is investigated using PHIL testing of the FESS and the real-time simulation of the studied microgrid. Both simulation and PHIL results confirm that the proposed adaptive inertia emulation controller can significantly reduce the maximum ROCOF and improve the frequency nadir during large frequency disturbances. Moreover, it outperforms several previously reported interval-based controllers in the investigated scenarios, while not requiring significantly more energy from the FESS.

# **8 Conclusions and Outlook**

# **8.1 Conclusions**

This work investigated the use of a FESS for the application of power smoothing and frequency support in low voltage distribution grids from various perspectives, including its allocation, sizing, modeling and model validation, real-time simulation, and PHIL testing in several scenarios.

Firstly, in order to identify the most suitable location for installing a FESS for power smoothing applications, a new data-driven method for estimating the relative voltage sensitivity to active power changes was introduced. The proposed method uses the mutual information score between the voltage and active power measurements to reflect the relative voltage sensitivity to active power changes at different locations of the grid, without the need of having a grid model. This method has been applied for allocating a FESS in a local grid in southern Germany with more than 1200 installation candidates. It was shown that the proposed method can successfully identify the locations, where the voltage is more sensitive to active power changes, and therefore, can benefit the most from a smoother power profile by using a FESS.

Next, a novel methodology for sizing energy storage systems based on historical measurement data was developed. The main contribution of the proposed method is the use of the motif discovery algorithm for detecting reoccurring consumption patterns in power profiles. The detected patterns are then used as the representative of the whole data set for deriving the energy storage characteristics. The method has been employed for sizing two different types of energy storage system, including a FESS for power smoothing applications. It has been demonstrated that the storage systems with the characteristics derived from only the detected patterns can fulfil their installation purposes such as power smoothing and peak shaving for most of the days throughout the whole the measurement period.

In the following chapters, a dynamic model for a high-speed FESS has been presented and validated with experimental results. The advantage of the presented model is that it incorporates the losses and the auxiliary power requirements of the FESS, which leads to an accurate estimation of the state of charge of the FESS during idle, charging, and discharging operational modes of the FESS, in all the investigated scenarios. The maximum difference between the state of charge of the developed model and the real FESS, among all investigated scenarios, has been approximately 0.8 % only.

For verifying the performance of the FESS, a laboratory setup for PHIL testing of a 60 kW high-speed FESS has been built and successfully tested. The setup allows testing the FESS in various scenarios and grid configurations for different applications, and also controlling and monitoring the FESS during the experiments. In order to evaluate the response of the FESS for the application frequency support using this setup, several frequency deviation scenarios including a grid code compliance verification test and the major frequency disturbance of August 9, 2019, in the power system of the UK have been simulated in real-time. It has been demonstrated that the FESS can respond very quickly to frequency deviations, and reaches the required power indicated by the grid code in just under 60 ms.

In order to demonstrate the advantages of the rapid response of the FESS for low-inertia power systems, the use of a FESS for providing emulated inertia in a low voltage microgrid has been investigated. A novel adaptive controller for inertia emulation using a FESS has been introduced, and its performance has been verified using simulations and experiments in two different scenarios in a low-inertia microgrid. The proposed controller combines the advantages of the fast reaction time of bang-bang controllers with the flexibility of self-adaptive ones in order to improve the inertia support provided by the FESS. It also considers the state of charge of the FESS and its power limitations at low rotational speeds into the control structure. Simulation results indicate that during the islanding of the studied microgrid, the use of the FESS with the proposed adaptive inertia emulation controller reduces the maximum ROCOF by 28 %, and the maximum frequency deviation by 44 %, in comparison to the case without the FESS. Also, the proposed controller has demonstrated a better performance in terms of reducing the maximum ROCOF and improving the frequency nadir in comparison to several previously reported adaptive controllers. The suggested controller has also been implemented on the 60 kW high-speed FESS using the concept of rapid control prototyping, and the performance of the FESS equipped with the proposed inertia emulation controller has been verified by means PHIL testing. A good match between the simulation and PHIL testing results has been observed, indicating the effectiveness of the proposed controller and the accuracy of the FESS model. However, slightly higher ROCOF values are observed during the PHIL testing, due to the minor delays present in the real FESS, which demonstrates the advantage of PHIL testing, where real component behaviours are reflected in the grid variables. It is important to note that the results presented in this work are the very first experimental validation of inertia emulation using a FESS.

In summary, this thesis provided a comprehensive study on the grid integration of a highspeed FESS, showing its advantages, such as a quick response time to power imbalances, and how new tools such as PHIL testing and data-driven methods can be used effectively for conducting grid integration studies.

# **8.2 Outlook**

There are several areas where the presented work in this thesis can be extended in the future. With the increasing number of measurement devices in power systems (e.g., smart meters, phasor measurement units), the data-driven approaches introduced for the allocation and sizing of FESS have a great potential to be applied to other use cases and applications. When measurement devices are available at different locations of the grid, the proposed allocation procedure based on the mutual information score can be fully automated, and applied for the allocation of other distributed energy resources. Furthermore, a new method has been recently introduced to estimate the mutual information score in real-time using a sliding window, by updating only a fixed number of nearest-neighbor relationships [175]. Therefore, the use of mutual information, as a measure of linear and nonlinear dependency between two variables, can be investigated for several other applications, such as real-time topology change detection (see [306]), improving load modelling, and adaptive controllers, such as the controller proposed in [169].

The use of the motif discovery algorithm for extracting reoccurring patterns in active power measurements can drastically reduce the input data size, and therefore, the computation time required for solving complex nonlinear optimization problems, which are commonly used for component sizing and system planning. This can be significantly beneficial, in particular when the optimization problem requires data over longer horizons, such as several years. As a future work, it would be interesting to quantify the impact of using motif discovery on the computation time of such optimization frameworks and on their outcome, and compare it to commonly used methods, such as clustering techniques.

The validated model of the FESS, which can accurately reflect its self-discharge rate, can be used by other researchers for simulating and investigating the advantages of a FESS in other applications, or designing and validation new controller. Since the model can be easily simulated in real-time using commercial real-time simulators, it can also be used for Controller Hardware-in-the-Loop (CHIL) testing of new controllers and energy management systems. Moreover, this model together with a power amplifier, that is capable to operate as a current source, can be used to build a FESS emulator, which is a cheaper and safer alternative of having the real FESS. The modeling can be improved in the future by disaggregating different losses and the power of each auxiliary component. This can help to build a more efficient FESS or operating it with less power requirements. Also, the small delays found within the real FESS can be added to the model to improve its accuracy in the millisecond range.

But more importantly, the laboratory setup established within this work for the PHIL testing of the FESS, can be used for evaluating the performance of a real FESS for other applications, including power smoothing and ramp rate control of renewables, voltage regulation, and supporting pulsed power loads in various low voltage grids and microgrids. In addition, the combination of the FESS with other storage technologies with a higher energy density, such as lithium-ion batteries, or with distributed energy resources which operate more efficiently under a constant operating point (e.g., fuel cells, micro gas turbines), can be investigated with the real FESS using this setup. As shown, this setup also enables to test new control strategies or energy management systems on the real FESS for each application, using the concept of rapid control prototyping. By combining rapid control prototyping and PHIL testing, this setup allows to verify the performance of these new controllers under the practical limitations and constraints of a real FESS, which are very often neglected in numerical simulation studies found in literature. These include auxiliary power requirements and losses of the FESS, its power limitations at low rotational speeds, delays and other system nonlinearities.

While the FESS has proven to be a viable solution for short-term high-power applications, further improvements in its efficiency (for instance, by using more efficient power converters, electrical machines, and bearings), reduction in its auxiliary power requirements, and a more competitive capital costs can potentially lead to an increased presence of this technology in power systems, in particular in low-inertia systems.

# **Bibliography**


and T. Leibfried, "Establishing a hardware-in-the-loop research environment with a hybrid energy storage system," in *IEEE PES Innovative Smart Grid Technologies Conference Europe*, Melbourne, 2016, pp. 497–503.


fluctuations in all-electric ship propulsion systems," *Applied Energy*, vol. 212, pp. 919–930, 2018.


2017.


utilizing a high-temperature superconducting bearing," *IEEE Transactions on Applied Superconductivity*, vol. 17, no. 2, pp. 2133–2137, 2007.


Chevrefils, M. Matar, R. Iravani, C. Dufour, J. Belanger, M. O. Faruque, K. Strunz, and J. A. Martinez, "Interfacing issues in real-time digital simulators," *IEEE Transactions on Power Delivery*, vol. 26, no. 2, pp. 1221–1230, 2011.


in-the-loop (PHIL) simulation," in *IEEE Power and Energy Society General Meeting*, National Harbor, MD, USA, 2014.


*Technologies Symposium (ESTS 2019)*, Washington, DC, USA, 2019, pp. 39–44.


Applications," *IEEE Transactions on Industrial Electronics*, vol. 65, no. 8, pp. 6612–6624, 2018.


*Pattern Analysis and Machine Intelligence*, vol. 27, no. 8, pp. 1226–1238, 2005.


incremental time series clustering," in *2013 IEEE International Conference on Big Data, Big Data 2013*, Silicon Valley, CA, USA, 2013, pp. 29–36.


*Industrial Electronics Conference (IECON)*, Florence, Italy, 2016.


Germany, 2014, pp. 153–160.


control of an IPM machine based on signal injection considering inductance saturation," *IEEE Transactions on Power Electronics*, vol. 28, no. 1, pp. 488–497, 2013.


extended Ward equivalent approach for power system security assessment," *Electric Power Systems Research*, vol. 42, no. 3, pp. 181–188, 1997.


*on Power Electronics*, vol. 33, no. 10, pp. 8488–8499, 2018.


8992–9002, 2019.


[307] IEC, "IEC 61000-2-2: Electromagnetic compatibility (EMC) - Part 2-2: Environment - Compatibility levels for low-frequency conducted disturbances and signalling in public low-voltage power supply systems," 2002.

# **Publications**

#### **Main publications**:


#### **Other Contributions:**

9. V. Arzamasov, R. Schwerdt, **S. Karrari**, K. Böhm, T. B. Nguyen, "Privacy Measures and Storage Technologies for Battery-Based Load Hiding—an Overview and Experimental Study," in Proceedings of the 11th ACM International Conference on Future Energy Systems (e-Energy 2020), Melbourne, 2020.


#### **Conference Poster and Oral Presentations:**


# **Appendix A. A Summary of Power Quality Measurements**

As described in subsection 4.1.4.2, in this work, measurement data have been collected from several low voltage distribution grids in southern Germany, which were used for the purpose of allocation and sizing of the FESS using data-driven methods. These data were recorded using A. Eberle's "PQI-DA smart" device, which is a power quality and disturbance recorder. In this appendix, after giving a short summary of the power quality requirements for the low voltage distribution grid in Germany, the status of the power quality at the measured low voltage grids are summarized.

In short, at none of the measured locations, a major power quality violation was observed, and a high power quality is being delivered by the distribution network operator at the investigated locations. However, several interesting incidents and observations were recorded, which are presented in this appendix.

# **A.1 Power Quality Requirements for Low Voltage Grids in Germany**

This section summarizes the minimum power quality requirements that the system operator should provide for the private and industrial customers, connected to the low voltage networks. These requirements include the maximum variation in the voltage amplitude, frequency, harmonics, flicker, and control signal level. These requirements are provided in the grid codes DIN EN 50160 and DIN EN 61000-2-2 for public grids (or DIN EN 61000-2-4 for industrial grids).

The objective of the system operator is to provide a balanced three-phase sinusoidal voltage with 50 Hz frequency. However, the variations in load, generation, switching of devices and lines, nonlinear loads, transformers during energizing, and faults may lead to deviation of these properties from their ideal values. Nevertheless, grid codes obligate the network operators to keep the electrical characteristics of the provided electricity supply within some predefined limits. Table A.1 summarizes the grid code requirements of the electricity supply in the low voltage networks in Germany. More details about each specific requirement can be found in DIN EN 50160.


Table A.1: Summary of the power quality requirements at the low voltage level in Germany [159], [307].

Figure A.1: Permissible relative voltage amplitudes of individual harmonics according the grid code DIN N 50160 [159].

Among the nine different low voltage distribution grids, at which measurements were conducted, no major power quality problem according to the criteria given in Table A.1 has been observed. As an example, the power quality indices measured at one of the low voltage distribution grids are presented in Table A.2.


Table A.2: Measured power quality indices at one of the low voltage distribution grids.

## **A.2 Summary of Interesting Recorded Incidents**

In this section, some interesting observations from the recorded measurement data are presented, regarding voltage deviations, asymmetries, and voltage and current harmonics.

In terms of the voltage magnitude, a voltage deviation larger than ±10 % of the nominal value has been barely witnessed during the measurement period. A rare example of a voltage drop of more than 13 % is shown in Figure A.2. However, such voltage drops were far from being a power quality issue, as they were very rare and lasted for a very short period.

Figure A.2: An example of voltage drop above ±10% of the nominal value. Such voltage drops were extremely rare in the measurement and often for a noticeably short time.

However, voltage deviations above the ±3 % limit, from which the distributed energy resources and energy storage systems should support the grid according to the grid code VDE-AR-N 4105:2018-11 (see Figure 4.5(b)), has been recorded quite often. This has been especially the case in low voltage distribution grid with a high share PV generation system, where overvoltage incidents has been mostly observed. Therefore, it can be concluded that the reactive power support by the distributed energy resources is currently required often within each day at the measured locations. With the increasing share of PV in these areas, more voltage deviations in these areas are expected. An example of a recorded over-voltage incident above 3 % at one of the low voltage distribution grids is shown in Figure A.3.

Figure A.3: An example of voltage deviations above ±3%, from which the distributed energy resources and energy storage system should support the grid with reactive power, according the German grid code VDE-ARN 4105:2018-11.

In terms of voltage and current asymmetry, voltage asymmetries were far below the power quality violation limit of 2 %, as discussed, reaching 0.67 % in the worst case at the measured locations. However, extremely high current asymmetries were recorded in almost all the measured low voltage grids, which has led to a high neutral current and negative sequence component. In some cases, the neutral current has been much greater than each individual phase currents, as illustrated in Figure A.4, as an example. Here, the negative sequence current has been multiple times more than the positive sequence current with the increase in the PV generation power around noon. This could be caused by old single-phase PV generation units, which increase the asymmetry in the grid around

noon. It should be noted that grid connection of new single-phase PV system is no longer permitted in Germany since several years.

Figure A.4: An example of high current asymmetry during one day in one of the measured low voltage distribution grids. (a) Relative negative sequence current, (b) Phases and neutral currents.

Regarding the voltage and current harmonics, an interesting and yet expected observation is the increase in the THD of the measured voltages with the increase in the PV generation power around noon, at the low voltage grid with a high share of PV generation. This can be clearly seen in the example illustrated in Figure A.5. Here, the increase in demand during the day coincides with the increase the PV generation, leading to a reduced current drawn from the medium voltage grid and the local consumption of the PV generation around noon. Therefore, the current drawn from the medium voltage grid is reduced at this time. However, the voltage THD continuous to increase with the increase in the PV generation. While the this is still far away from the 8 % threshold to be a power quality concern, with the increase in the PV installations at these low voltage grids, harmonics can become an issue in the future. An example of currents drawn at such locations is shown in Figure A.6.

Figure A.5: An example of the increase in the THD of the voltage at noon with increase in the PV generation.

Figure A.6: An example of current and voltage waveforms at a low voltage distribution grid with a high THD.

# **Appendix B. Parameters of the Microgrid Components**

The parameters of the distributed energy resources added to the CIGRE European low voltage distribution grid benchmark in Chapter 7 are presented in Tables B.1 and B.2.


Table B.1: Parameters of the diesel-based Combined Heat and Power (CHP) system.


Table B.2: Parameters of the photovoltaic systems.

# **Table of Figures**









# **Table of Tables**




# **Abbreviations**


With the aim of decarbonizing the electricity sector, the share of renewable energy sources in power systems around the globe is consistently increasing. However, due to the intermittent nature of these sources, maintaining the instantaneous balance between the generation and the demand, and therefore, the grid frequency, can be challenging. In addition, since the converter-interfaced renewables do not inherently provide inertia to the system, the cumulative system inertia is simultaneously decreasing, resulting in faster changes in the grid frequency. A Flywheel Energy Storage System (FESS) can rapidly inject or absorb high amounts of active power in order to support the grid, following an abrupt change in the generation or in the demand. In addition to the quick response time, a FESS has the advantage of a high power density and a large number of cycles, with no capacity loss throughout its lifetime. These characteristics make the FESS a well-suited candidate for providing frequency support to the grid or smoothing out short-term power variations at a local level.

The work presented in this book studies the grid integration of a high-speed FESS in low voltage distribution grids from several perspectives, including optimal allocation, sizing, modeling, real-time simulation, and Power Hardware-in-the-Loop (PHIL) testing. Moreover, this work presents the very first experimental validation of inertia emulation using a FESS.

INTEGRATION OF FLYWHEEL ENERGY STORAGE SYSTEMS IN LOW VOLTAGE DISTRIBUTION GRIDS

Shahab Karrari